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THE INTERPOLATION THEORY 
OF RADIAL BASIS FUNCTIONS 

B. J. C. Baxter 
Summary 

The problem of interpolating functions of d real variables (d > 1) occurs naturally 
in many areas of applied mathematics and the sciences. Radial basis function 
methods can provide interpolants to function values given at irregularly positioned 
points for any value of d. Further, these interpolants are often excellent approxi- 
mations to the underlying function, even when the number of interpolation points 
is small. 

In this dissertation we begin with the existence theory of radial basis function 
interpolants. It is first shown that, when the radial basis function is a p-norm 
and 1 < p < 2, interpolation is always possible when the points are all different 
and there are at least two of them. Our approach extends the analysis of the case 
p = 2 devised in the 1930s by Schoenberg. We then show that interpolation is not 
always possible when p > 2. Specifically, for every p > 2, we construct a set of 
different points in some lZ d for which the interpolation matrix is singular. This 
construction seems to have no precursor in the literature. 

The greater part of this work investigates the sensitivity of radial basis func- 
tion interpolants to changes in the function values at the interpolation points. 
This study was motivated by the observation that large condition numbers occur 
in some practical calculations. Our early results show that it is possible to recast 
the work of Ball, Narcowich and Ward in the language of distributional Fourier 
transforms in an elegant way. We then use this language to study the interpola- 
tion matrices generated by subsets of regular grids. In particular, we are able to 
extend the classical theory of Toeplitz operators to calculate sharp bounds on the 
spectra of such matrices. Moreover, we also describe some joint work with Charles 
Micchelli in which we use the theory of Polya frequency functions to continue this 
work, as well as shedding new light on some of our earlier results. 
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Applying our understanding of these spectra, we construct preconditioners for 
the conjugate gradient solution of the interpolation equations. The preconditioned 
conjugate gradient algorithm was first suggested for this problem by Dyn, Levin 
and Rippa in 1983, who were motivated by the variational theory of the thin 
plate spline. In contrast, our approach is intimately connected to the theory of 
Toeplitz forms. Our main result is that the number of steps required to achieve 
solution of the linear system to within a required tolerance can be independent of 
the number of interpolation points. In other words, the number of floating point 
operations needed for a regular grid is proportional to the cost of a matrix-vector 
multiplication. The Toeplitz structure allows us to use fast Fourier transform 
techniques, which implies that the total number of operations is a multiple of 
nlogn, where n is the number of interpolation points. 

Finally, we use some of our methods to study the behaviour of the multiquadric 
when its shape parameter increases to infinity. We find a surprising link with 
the sinus cardinalis or sine function of Whittaker. Consequently, it can be highly 
useful to use a large shape parameter when approximating band-limited functions. 
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1 : Introduction 

The multivariate interpolation problem occurs frequently in many branches of 
science and engineering. Typically, we are given a discrete set / in TZ d , where d is 
greater than one, and real numbers {f%}i£i- Our task is to construct a continuous 
or sufficiently differentiable function s: 1Z d — > 1Z such that 

s(i) = fi, i E I, (1.1) 

and we say that s interpolates the data {(z, fi):ie I}. Interpolants can be highly 
useful. For example, we may need to approximate a function whose values are 
known only at the interpolation points, that is we are ignorant of its behaviour 
outside /. Alternatively, the underlying function might be far too expensive to 
evaluate at a large number of points, in which case the aim is to choose an in- 
terpolant which is cheap to compute. We can then use our interpolant in other 
algorithms in order to, for example, calculate approximations to extremal values of 
the original function. Another application is data-compression, where the size of 
our initial data {(i, fi):i<E 1} exceeds the storage capacity of available computer 
hardware. In this case, we can choose a subset I of I and use the corresponding 
data to construct an interpolant with which we estimate the remaining values. It 
is important to note that in general I will consist of scattered points, that is its 
elements can be irregularly positioned. Thus algorithms that apply to arbitrary 
distributions of points are necessary. Such algorithms exist and are well under- 
stood in the univariate case (see, for instance, Powell (1981)), but many difficulties 
intrude when d is bigger than one. 

There are many applications of multivariate interpolation, but we prefer to 
treat a particular application in some detail rather than provide a list. Therefore 
we consider the following interesting example of Barrodale et al (1991). 

When a time-dependent system is under observation, it is often necessary 
to relate pictures of the system taken at different times. For example, when 
measuring the growth of a tumour in a patient, we must expect many changes to 
occur between successive X-ray photographs, such as the position of the patient 
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or the amount of fluid in the body's tissues. If we can identify corresponding 
points on the two photographs, such as parts of the bone structure or intersections 
of particular veins, then these pairs of points can be viewed as the data for two 
interpolation problems. Specifically, let (xj, yj)?=i be the coordinates of the points 
in one picture, and let the corresponding points in the second picture be ?7j)" =1 . 
We need functions s x :1Z 2 — >■ 1Z and s y :1Z 2 — >■ 1Z such that 

s x (xj, Vj) = €j and s v (xj, yj) = rjj for j = 1, . . . , n. (1.2) 

Therefore we see that the scattered data interpolation problem arises quite nat- 
urally as an attempt to approximate the non-linear coordinate transformation 
mapping one picture into the next. 

It is important to understand that interpolation is not always desirable. For 
example, our data may be corrupted by measurement errors, in which case there 
is no good reason to choose an approximation which satisfies the interpolation 
equations, but we do want to construct an approximation which is close to the 
function values in some sense. One option is to choose our function s:lZ d — > 71 
from some family (usually a linear space) of functions so as to minimize a certain 
functional G, such as 

G(s-f) = Y J [h-sm 2 , (1-3) 

iei 

which is the familiar least-squares fitting problem. Of course this can require 
the solution of a nonlinearly constrained optimization problem, depending on the 
family of functions and the functional G. Another alternative to interpolation 
takes s to be the sum of decaying functions, each centred at a point in / and 
taking the function value at that point. Such an approximation is usually called 
a quasi-interpolant, reflecting the requirement that it should resemble the inter- 
polant in some suitable way. These methods are of both practical and theoretical 
importance, but we emphasize that this dissertation is restricted to interpolation, 
specifically interpolation using radial basis functions, for which we refer the reader 
to Section 1.5 and the later chapters of the dissertation. 
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We now briefly describe some other multivariate approximation schemes. Of 
course, our treatment does not provide a thorough overview of the field, for which 
we refer the reader to de Boor (1987), Franke (1987) or Hayes (1987). However, 
it is interesting to contrast radial basis functions with some of the other methods. 
In fact, the memoir of Franke (1982) is dedicated to this purpose; it contains 
careful numerical experiments using some thirty methods, including radial basis 
functions, and provides an excellent reason for their theoretical study: they obtain 
excellent accuracy when interpolating scattered data. Indeed, Franke found them 
to excel in this sense when compared to the other tested methods, thus providing 
an excellent reason for their theoretical study. 

1.1 Polynomial interpolation 

Let P be a linear space of polynomials in d real variables spanned by {p%)i^i, where 
I is the discrete subset of 1Z d discussed at the beginning of the introduction. Then 
an interpolant s: 1Z d — >■ 1Z of the form 



exists if and only if the matrix (pi(j))i,jei is invertible. We see that this prop- 
erty depends on the geometry of the centres when d > 1, which is a signifi- 
cant difficulty. One solution is to choose a particular geometry. As an example 
we describe the tensor product approach on a "tartan grid". Specifically, let 
/ = {(xj, yk) '■ 1 < j < 1,1 < k < m}, where x\ < ■ • • < x\ and yi < ■ • • < y m 
are given real numbers, and let {f( x -,y k ) '■ 1 < j < 1,1 < k < m} be the func- 
tion values at these centres. We let (Lj)^ =1 and {L\)™ =1 be the usual univariate 
Lagrange interpolating polynomials associated with the numbers (xj)\ and {yk)T 
respectively and define our interpolant s: 1Z 2 — > 1Z by the equation 

l m 

S (x, y) = J2J2 k^,vu) L )i x ) L l{yl y) e n 2 . (1.5) 
j=i k=i 

Clearly this approach extends to any number of dimensions d. 





(1.4) 
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1.2 Tensor product methods 

The tensor product scheme for tartan grids described in the previous section is 
not restricted to polynomials. Using the same notation as before, we replace 
(Lj) l j =1 and (L\)™ =1 by univariate functions (Pj) l j =1 and {Qk)k=i respectively. 
Our interpolant takes the form 

l m 

s(x,y) = ^2^2y jk Pj(x)Q k (y), (x,y)eTZ 2 , (1.7) 

j=l k=l 

from which we obtain the coefficients (yjk)- By adding points outside the interval 
[xi, xi] and [yi,y m ] we can choose (Pj) and (Qk) to be univariate B-splines. In this 
case the linear systems involved are invertible and banded, so that the number of 
operations and the storage required are both multiples of the total number of points 
in the tartan grid. Such methods are extremely important for the subtabulation 
of functions on regular grids, and clearly the scheme exists for any number of 
dimensions d. A useful survey is the book of Light and Cheney (1986) 

1.3 Multivariate Splines 

Generalizing some of the properties of univariate splines to a multivariate set- 
ting has been an idee fixe of approximation theory. Thus the name "spline" is 
overused, being applied to almost any extension of univariate spline theory. In this 
section we briefly consider box splines. These are compactly supported piecewise 
polynomial functions which extend Schoenberg's characterization of the S-spline 
B{-; to, . . . , tk) with arbitrary knots to, . . . , tk as the "shadow" of a /c-dimensional 
simplex (Schoenberg (1973), Theorem 1, Lecture 1). Specifically, the box spline 
B(-;A) associated with the d x n matrix A is the distibution defined by 

£(•; A) : C^(TZ d ) ->■ K : cp ^ / tp(Ax) dx, 

where C Q yD (1Z d ) is the vector subspace of C°°{1Z d ) whose elements vanish at infinity. 
If we let a\, . . . , a n e 1Z d be the columns of A, then the Fourier transform of the 
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box spline is given by 

n 

where sinc(x) = sin(x/2)/(x/2). We see that a simple example of a box spline 
is a tensor product of univariate B-splines. It can be shown that there exist box 
splines with smaller supports than tensor product B-splines. 

A large body of mathematics now exists, and a suitable comprehensive 
review is the long paper of Dahmen and Micchelli (1983). Further, this theory is 
also yielding useful results in the study of wavelets (see Chui (1992)). However, 
there are many computational difficulties. At present, box spline software is not 
available from the main providers of scientific computation packages. 

1.4 Finite element methods 

Finite element methods can provide extremely flexible piecewise polynomial spaces 
for approximation and scattered data interpolation. When d = 2 we first choose 
a triangulation of the points. Then a polynomial is constructed on each triangle, 
possibly using function values and partial derivative values at other points in 
addition to the vertices of the triangulation. This is a non-trivial problem, since 
we usually require some global differentiability properties, that is the polynomials 
must fit together in a suitably smooth way. Further, the partial derivatives are 
frequently unknown, and these methods can be highly sensitive to the accuracy of 
their estimates (Franke (1982)). 

Much recent research has been directed towards the choice of triangula- 
tion. The Delaunay triangulation (Lawson (1977)) is often recommended, but 
some work of Dyn, Levin and Rippa (1986) indicates that greater accuracy can be 
achieved using data- dependent triangulations, that is triangulations whose com- 
ponent triangles reflect the geometry of the function in some way. Finally, the 
complexity of constructing triangulations in higher dimensions effectively limits 
these methods to two and three dimensional problems. 
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1.5 Radial basis functions 

A radial basis function approximation takes the form 

s (x) = J2yM\\x-i\\), xen d , (1.8) 

iei 

where ip: [0, oo) — > 1Z is a fixed univariate function and the coefficients (yi)iei 
are real numbers. We do not place any restriction on the norm || ■ || at this 
point, although we note that the Euclidean norm is the most common choice. 
Therefore our approximation s is a linear combination of translates of a fixed 
function x i-> ^(||x||) which is "radially symmetric" with respect to the given 
norm, in the sense that it clearly possesses the symmetries of the unit ball. We 
shall often say that the points {xj)™ =1 are the centres of the radial basis function 
interpolant. Moreover, it is usual to refer to <p> as the radial basis function, if the 
norm is understood. 

If / is a finite set, say I = (xj)" =1 , the interpolation conditions provide the 
linear system 

Ay = /, (1.9) 

where 

A=U\\x J -x k \\)) n (1.10) 

y = ( yj )? =1 and / = (f 3 )] =1 . 

One of the most attractive features of radial basis function methods is the 
fact that a unique interpolant is often guaranteed under rather mild conditions 
on the centres. In several important cases, the only restrictions are that there 
are at least two centres and they are all distinct, which are as simple as one 
could wish. However, one important exception to this statement is the thin plate 
spline introduced by Duchon (1975, 1976), where we choose (p(r) = r 2 logr. It 
is easy to see that the interpolation matrix A given by (1.10) can be singular for 
non-trivial sets of distinct centres. For example, choosing a^2, • • • , x n to be any 
different points on the sphere of unit radius whose centre is x\, we conclude that 
the first row and column of A consist entirely of zeros. Of course, such examples 
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exist for any function <p with more than one zero. Fortunately, it can be shown 
that it is suitable to add a polynomial of degree m > 1 to the definition of s 
if the centres are unisolvent, which means that the zero polynomial is the only 
polynomial of degree m which vanishes at every centre (see, for instance, Powell 
(1992)). The extra degrees of freedom are usually taken up by moment conditions 
on the coefficients (yj)™ =1 . Specifically, we have the equations 

n 

^2y k ip(\\xj -x k \\) + P(xj) = fj, j = 1,2, 

k=1 (1.11) 

^VkP{xk) = for every p G U m (n d ), 
k=i 

where Um(TZ d ) denotes the vector space of polynomials in d real variables of total 
degree m, and the theory guarantees the existence of a unique vector (y 3 -)?=i an d 
a unique polynomial P G Ii m {lZ d ) satisfying (1.11). Moreover, because (1.8) does 
not reproduce polynomials when / is a finite set, it is sometimes useful to augment 
s in this way. 

In fact Duchon derived (1.11) as the solution to a variational problem when 
d = 2: he proved that the function s given by (1.11) minimizes the integral 

Jtz 2 

where m = 1 and s satisfies some differentiability conditions. Duchon's treatment 
is somewhat abstract, using sophisticated distribution theory techniques, but a 
detailed alternative may be found in Powell (1992). We do not study the thin 
plate spline in this dissertation, although many of our results are highly relevant 
to its behaviour. 

In his comparison of multivariate approximation methods, Franke (1982) 
considered several radial basis functions including the thin plate spline. Therefore 
we briefly consider some of these functions. 

The multiquadric 

Here we choose ip(r) = (r 2 + c 2 ) 1 / 2 , where c is a real constant. The interpolation 
matrix A is invertible provided only that the points are all different and there are 

7 



Introduction 

at least two of them. Further, this matrix has an important spectral property: it 
is almost negative definite; we refer the reader to Section 2 for details. 

Franke found that this radial basis function provided the most accurate 
interpolation surfaces of all the methods tried for interpolation in two dimensions. 
His centres were mildly irregular in the sense that the range of distances between 
centres was not so large that the average distance became useless. He found that 
the method worked best when c was chosen to be close to this average distance. 
It is still true to say that we do not know how to choose c for a general function. 
Buhmann and Dyn (1991) derived error estimates which indicated that a large 
value of c should provide excellent accuracy. This was borne out by some calcu- 
lations and an analysis of Powell (1991) in the case when the centres formed a 
regular grid in one dimension. Specifically, he found that the uniform norm of the 
error in interpolating f(x) = x 2 on the integer grid decreased by a factor of 10 3 
when c increased by one; see Table 6 of Powell (1991) for these stunning results. 
In Chapter 7 of this thesis we are able to show that the interpolants converge uni- 
formly as c — > oo if the underlying function is square-integrable and band-limited, 
that is its Fourier transform is supported by the interval [— tt, iv] d . Thus, for many 
functions, it would seem to be useful to choose a large value of c. Unfortunately, 
if the centres form a finite regular grid, then we find that the smallest eigenvalue 
of the distance decreases exponentially to zero as c tends to infinity. Indeed, the 
reader is encouraged to consider Table 4.1, where we find that the smallest eigen- 
value decreases by a factor of about 20 when c is increased by one and the spacing 
of the regular grid is unity. 

We do not consider the polynomial reproduction properties of the multi- 
quadric discovered by Buhmann (1990) in this dissertation, but we do make use 
of some of his work, in particular his formula for the cardinal function's Fourier 
transform (see Chapter 7). However, we cannot resist mentioning one of the bril- 
liant results of Buhmann, in particular the beautiful and surprising result that 
the degree of polynomials reproduced by interpolation on an infinite regular grid 
actually increases with the dimension. The work of Jackson (1988) is also highly 
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relevant here. 
The Gaussian 

There are many reasons to advise users to avoid the Gaussian (p(r) = exp(— cr 2 ). 
Franke (1982) found that it is very sensitive to the choice of parameter c, as we 
might expect. Further, it cannot even reproduce constants when interpolating 
function values given on an infinite regular grid (see Buhmann (1990)). Thus 
its potential for practical computer calculations seems to be small. However, 
it possesses many properties which continue to win admirers in spite of these 
problems. In particular, it seems that users are seduced by its smoothness and 
rapid decay. Moreover the Gaussian interpolation matrix (1.10) is positive definite 
if the centres are distinct, as well as being suited to iterative techniques. I suspect 
that this state of affairs will continue until good software is made available for 
radial basis functions such as the multiquadric. Therefore I wish to emphasize 
that this thesis addresses some properties of the Gaussian because of its theoretical 
importance rather than for any use in applications. 

In a sense it is true to say that the Gaussian generates all of the radial 
basis functions considered in this thesis. Here we are thinking of the Schoenberg 
characterization theorems for conditionally negative definite functions of order 
zero and order one. These theorems and related results occur many times in this 
dissertation. 

The inverse multiquadric 

Here we choose (p(r) = (r 2 + c 2 ) -1 / 2 . Again , Franke (1982) found that this radial 
basis function can provide excellent approximations, even when the number of 
centres is small. As for the multiquadric, there is no good choice of c known at 
present. However, the work presented in Chapter 7 does extend to this function 
(although this analysis is not presented here), so that sometimes a large value of 
c can be useful. 

The thin plate spline 

We have hardly touched on this highly important function, even though the works 

9 



Introduction 

of Franke (1982) and Buhmann (1990) indicate its importance is two dimensions 
(and, more generally, in even dimensional spaces). However, we aim to generalize 
the norm estimate material of Chapters 3-5 to this function in future. There 
is no numerical evidence to indicate that this ambition is unfounded, and the 
preconditioning technique of Chapter 6 works equally well when applied to this 
function. Therefore we are optimistic that these properties will be understood 
more thoroughly in the near future. 

1.6 Contents of the thesis 

Like Gaul, this thesis falls naturally into three parts, namely Chapter 2, Chapters 
3-6, and Chapter 7. In Chapter 2 we study and extend the work of Schoenberg 
and Micchelli on the nonsingularity of interpolation matrices. One of our main 
discoveries is that it is sometimes possible to prove nonsingularity when the norm is 
non-Euclidean. Specifically, we prove that the interpolation matrix is non-singular 
if we choose a p-norm for 1 < p < 2 and if the centres are different and there are 
at least two of them. This complements the work of Dyn, Light and Cheney 
(1991) which investigates the case when p = 1. They find that a necessary and 
sufficient condition for nonsingularity when d = 2 is that the points should not 
form the vertices of a closed path, which is a closed polygonal curve consisting of 
alternately horizontal and vertical arcs. For example, the 1-norm interpolation 
matrix generated by the vertices of any rectangle is singular. Therefore it may 
be useful that we can avoid these difficulties by using a p-norm for some p G 
(1,2). However, the situation is rather different when p > 2. This is probably 
the most original contribution of this section, since it makes use of a device that 
seems to have no precursor in the literature and is wholly independent of the 
Schoenberg-Micchelli corpus. We find that, if both p and the dimension d exceed 
two, then it is possible to construct sets of distinct points which generate a singular 
interpolation matrix. It is interesting to relate that these sets were suggested by 
numerical experiment, and the author is grateful to M. J. D. Powell for the use of 
his TOLMIN optimization software. 
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The second part of this dissertation is dedicated to the study of the spectra 
of interpolation matrices. Thus, having studied the nonsingularity (or otherwise) 
of certain interpolation matrices, we begin to quantify . This study was initiated by 
the beautiful papers of Ball (1989), and Narcowich and Ward (1990, 1991), which 
provided some spectral bounds for several functions, including the multiquadric. 
Our main findings are that it is possible to use Fourier transform methods to 
address these questions, and that, if the centres form a subset of a regular grid, 
then it is possible to provide a sharp upper bound on the norm of the inverse of the 
interpolation matrix. Further, we are able to understand the distribution of all the 
eigenvalues using some work of Grenander and Szego (1984). This work comprises 
Chapters 3 and 4. In the latter section, it turns out that everything depends on an 
infinite product expansion for a Theta function of Jacobi type. This connection 
with classical complex analysis still excites the author, and this excitement was 
shared by Charles Micchelli. Our collaboration, which forms Chapter 5, explores 
a property of Polya frequency functions which generalizes the product formula 
mentioned above. Furthermore, Chapter 5 contains several results which attack 
the norm estimate problem of Chapter 4 using a slightly different approach. We 
find that we can remove some of the assumptions required at the expense of a 
little more abstraction. This work is still in progress, and we cannot yet say 
anything about the approximation properties of our suggested class of functions. 
We have included this work because we think it is interesting and, perhaps more 
importantly, new mathematics is frequently open-ended. 

Chapters 6 and 7 apply the work of previous chapters. In Chapter 6 we use 
our study of Toeplitz forms in Chapter 4 to suggest a preconditioner for the conju- 
gate gradient solution of the interpolation equations, and the results are excellent, 
although they only apply to finite regular grids. Of course it is our hope to extend 
this work to arbitrary point sets in future. We remark that our approach is rather 
different from the variational heuristic of Dyn, Levin and Rippa (1986), which 
concentrated on preconditioners for thin plate splines in two dimensions. Proba- 
bly our most important practical finding is that the number of iterations required 

11 



Introduction 

to attain a solution to within a particular tolerance seems to be independent of 
the number of centres. 

Next, Chapter 7 is unique in that it is the only chapter of this thesis 
which concerns itself with the approximation power of radial basis function spaces. 
Specifically, we investigate the behaviour of interpolation on an infinite regular grid 
using a multiquadric (p(r) = (r 2 + c 2 ) 1 / 2 when the parameter c tends to infinity. 
We find an interesting connection with the classical theory of the Whittaker car- 
dinal spline: the Fourier transform of the cardinal (or fundamental) function of 
interpolation converges (in the L 2 norm) to the characteristic function of the cube 
[— n,n] d . This enables us to show that the interpolants to certain band-limited 
functions converge uniformly to the underlying function when c tends to infinity. 

An aside Finally, we cannot resist the following excursion into the theory of 
conic sections, whose only purpose is to lure the casual reader. Let S and S' be 
different points in 1Z 2 and let f:1Z 2 — > TZ be the function defined by 

f(x) = \\x-S\\ + \\x-S% xen 2 , 

where || ■ || is the Euclidean norm. Thus the contours of / constitute the set of 
all ellipses whose focal points are S and S'. By direct calculation we obtain the 
expression 

™-(iHi) + (iKfi) 

which implies the relations 

whose geometric interpretation is the reflector property of the ellipse. A similar 
derivation exists for the hyperbola. 

1.7 Notation 

We have tried to use standard notation throughout this thesis with a few excep- 
tions. Usually we denote a finite sequence of points in <i-dimensional real space 
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lZ d by subscripted variables, for example (x_j)" =1 . However we have avoided this 
usage when coordinates of points occur. Thus Chapters 2 and 5 use superscripted 
variables, such as (x J )™ =1 , and coordinates are then indicated by subscripts. For 
example, x J k denotes the fcth coordinate of the jth vector of a sequence of vectors 
(a>0j=i- The inner product of two vectors x and y is denoted xy in the context 
of a Fourier transform, but we have used the more traditional linear algebra form 
x T y in Chapter 6 and in a few other places. We have used no special notation for 
vectors, and we hope that no ambiguity arises thereby. 

Given any absolutely integrable function /: lZ d — > 1Z, we define its Fourier 
transform by the equation 



We also use this normalization when discussing distributional Fourier transforms. 
Thus, if it is permissible to invert the Fourier transform, then the integral takes 
the form 



The norm symbol (|| ■ ||) will usually denote the Euclidean norm, but this 
is not so in Chapter 1. Here the Euclidean norm is denoted by | • | to distinguish 
it from other norm symbols. 

Finally, the reader will find that the term "radial basis function" can often 
mean the univariate function <p: [0, oo) — > 1Z and the multivariate function lZ d 3 
x !->■ <p(||x||). This abuse of notation was inherited from the literature and seems 
to have become quite standard. However, such potential for ambiguity is bad. It 
is perhaps unusual for the author of a dissertation to deride his own notation, but 
it is hoped that the reader will not perpetuate this terminology. 





x e n d . 
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2 : Conditionally positive functions and 
p-norm distance matrices 

2.1. Introduction 

The real multivariate interpolation problem is as follows. Given distinct points 
x 1 ,...,x n G 1Z d and real scalars /i,...,/ n , we wish to construct a continuous 
function s: 1Z d —¥ 1Z for which 

s{x l ) = fi, for % = 1, . . .,ra. 

The radial basis function approach is to choose a function </?: [0, oo) — > [0, oo) and 
a norm || • || on 1Z d and then let s take the form 

n 

s{x) = Y J \<P{h-x i \\)- 

i=l 

Thus s is chosen to be an element of the vector space spanned by the functions 
£ h-><£>(||£ — for i = 1, ... ,n. The interpolation conditions then define a linear 
system A A = /, where A G TZ nXn is given by 

^•=^(11^-^11), for 1 < < n, 

and where A = (Ai, A n ) and / = (/i, / n ). In this thesis, a matrix such as A 
will be called a distance matrix. 

Usually || ■ || is chosen to be the Euclidean norm, and in this case Micchelli 
(1986) has shown the distance matrix generated by distinct points to be invertible 
for several useful choices of <p. In this chapter, we investigate the invertibility of 
the distance matrix when || • || is a p-norm for 1 < p < oo, p ^ 2, and <p(t) = t, 
the identity. We find that p-norms do indeed provide invertible distance matrices 
given distinct points, for 1 < p < 2. Of course, p = 2 is the Euclidean case 
mentioned above and is not included here. Now Dyn, Light and Cheney (1991) 
have shown that the 1— norm distance matrix may be singular on quite innocuous 
sets of distinct points, so that it might be useful to approximate || • ||i by || • || p for 
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some p G (1,2]. This work comprises section 2.3. The framework of the proof is 
very much that of Micchelli (1986). 

For every p > 2, we find that distance matrices can be singular on certain 
sets of distinct points, which we construct. We find that the higher the dimension 
of the underlying vector space for the points x 1 , . . . , x n , the smaller the least p for 
which there exists a singular p-norm. 

2.2. Almost negative matrices 

Almost every matrix considered in this section will induce a non-positive form on 
a certain hyperplane in 1Z n . Accordingly, we first define this ubiquitous subspace 
and fix notation. 

Definition 2.2.1. For any positive integer n, let 

n 

Z n = {yeK n :J2vi = Q }- 

i=i 

Thus Z n is a hyperplane in 1Z n . We note that Z\ = {0}. 

Definition 2.2.2. We shall call A G 1Z nXn almost negative definite (AND) if A 
is symmetric and 

y 1 Ay < 0, whenever y G Z n . 

Furthermore, if this inequality is strict for all non-zero y G Z n , then we shall call 
A strictly AND. 

Proposition 2.2.3. Let A G 7Z nXn be strictly AND with non-negative trace. Then 

(-l) n_1 detA > 0. 

Proof. We remark that there are no strictly AND lxl matrices, and hence n > 2. 
Thus A is a symmetric matrix inducing a negative-definite form on a subspace of 
dimension n — 1 > 0, so that A has at least n — 1 negative eigenvalues. But trace 
A > 0, and the remaining eigenvalue must therefore be positive. | 
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Micchelli (1986) has shown that both A^ = \x* - x j \ and A^ = (1 + \x l - x j \ 2 )^ 
are AND, where here and subsequently | • | denotes the Euclidean norm. In fact, 
if the points distinct and n > 2, then these matrices are strictly 

AND. Thus the Euclidean and multiquadric interpolation matrices generated by 
distinct points satisfy the conditions for proposition 2.2.3. 

Much of the work of this chapter rests on the following characterization 
of AND matrices with all diagonal entries zero. This theorem is stated and used 
to good effect by Micchelli (1986), who omits much of the proof and refers us to 
Schoenberg (1935). Because of its extensive use we include a proof for the conve- 
nience of the reader. The derivation follows the same lines as that of Schoenberg 
(1935). 

Theorem 2.2.4. Let A G TZ nXn have all diagonal entries zero. Then A is AND 
if and only if there exist n vectors y 1 , . . . ,y n G lZ n for which 

A ij = \y i -yi\ 2 . 

Proof. Suppose A^ = \y l — y J '\ 2 for vectors y 1 , . . . , y n G 7Z n . Then A is symmetric 
and the following calculation completes the proof that A is AND. Given any z G 
Z n , we have 



Az= J2 z i z M -y J \ 2 

n 

= ^z i z j (\y i \ 2 + \yi\ 2 -2(y*) T (yi)) 



n 

—2 ^ ZiZj(y l ) T (yi) since the coordinates of z sum to zero, 

M = l 
n 

i=l 



2 

< 0. 



This part of the proof is given in Micchelli (1986). The converse requires two 
lemmata. 
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Lemma 2.2.5. Let B G TZ kxk be a symmetric non-negative definite matrix. Then 
we can find £ 1 , . . . , £ fc G lZ k such that 

B l3 = \e\ 2 + \e\ 2 -\e-t j \ 2 - 

Proof. Since B is symmetric and non-negative definite, we have B = P T P, for 
some P G 1Z kxk . Let p 1 , . . . , p fc be the columns of P. Thus 



Now 



Hence 



B ij = ( P r( P >). 



\p l -pJ\ 2 = \p l \ 2 + \p 3 \ 2 -2{p l ) T {p J ] 



All that remains is to define £ l = p 1 / ^/2 , for % = 1, . . . , k. | 

Lemma 2.2.6. Let A G 1Z nXn . Let e 1 , . . . , e n denote the standard basis for 1Z n , 
and define 

f % = e n -e\ fori = l,...,n—l, 
f n = e n . 

Finally, let F G 7?. nXn 6e £/ie matrix with columns f 1 ,...,f n . Then 

{-F T AF) lJ = A in + A nj - Aij - A nn , fori <i,j <n- 1, 

( P ^4P)in = 

(-P T AF) m = A m - A nn , for 1 < i < n - 1, 

( P ^4P)nn = 

Proo/. We simply calculate (-F T AF) y = -(f) T A(p). | 

We now return to the proof of Theorem 2.2.4: Let A G TZ nXn be AND with 
all diagonal entries zero. Lemma 2.2.6 provides a convenient basis from which 
to view the action of A. Indeed, if we set B = —F T AF, as in Lemma 2.2.6, 
we see that the principal submatrix of order n — 1 is non-negative definite, since 
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J 1 ,...,/" -1 form a basis for Z n . Now we appeal to Lemma 2.2.5, obtaining 
f 1 ,...,^" 1 G TZ 71 - 1 such that 

B tj = |f| 2 + l^l 2 - \C - , for 1 < ij < n - 1, 

while Lemma 2.2.6 gives 

Setting i = j and recalling that An = 0, we find 

A in = \C\ 2 , forl<z<n-l 

and thus we obtain 

Ar :i = \C-e' 2 - for l<i,j<n-l. 

Now define £ n = 0. Thus A i3 = |f - £ J | 2 , for 1 < i, j < n, where 
G TZ n ~ 1 . We may of course embed 7£ n_1 in 7£ n . More formally, 
let i-.VJ 1-1 <— > lZ n be the map t: (x\, . . . , £ n -i) >->■ (xi, . . . , x n _i, 0), and, for 
z = 1, . . . , n, define y l = t(C). Thus y 1 , . . . , y n G Tl n and 

A y = |y*-^| 2 . 

The proof is complete. | 

Of course, the fact that y n = by this construction is of no import; we 
may take any translate of the n vectors y , . . . , y n if we wish. 

2.3. Applications 

In this section we introduce a class of functions inducing AND matrices and then 
use our characterization Theorem 2.2.4 to prove a simple, but rather useful, the- 
orem on composition within this class. We illustrate these ideas in examples 
2.3.3-2.3.5. The remainder of the section then uses Theorems 2.2.4 and 2.3.2 to 
deduce results concerning powers of the Euclidean norm. This enables us to derive 
the promised p-norm result in Theorem 2.3.11. 
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Definition 2.3.1. We shall call /: [0, oo) — > [0, oo) a conditionally negative defi- 
nite function of order 1 ( CND1 ) if, for any positive integers n and d, and for any 
points x 1 , . . . , x n G TZ d , the matrix A G TZ nXn defined by 

Atj = f(\x l - x J \ 2 ), for 1 < i,j < n, 

is AND. Furthermore, we shall call / strictly CND1 if the matrix A is strictly 
AND whenever n > 2 and the points x 1 , . . . , x n are distinct. 

This terminology follows that of Micchelli (1986), Definition 2.3.1 . We see 
that the matrix A of the previous definition satisfies the conditions of proposition 
2.2.3 if / is strictly CND1, n > 2 and the points x 1 , . . . , x n are distinct. 

Theorem 2.3.2. 

(1) Suppose that f and g are CND1 functions and that /(0) = 0. Then go f is 
also a CND1 function. Indeed, if g is strictly CND1 and f vanishes only at 
0, then g o f is strictly CND1. 

(2) Let A be an AND matrix with all diagonal entries zero. Let g be a CND1 
function. Then the matrix defined by 

B ij = g( A ij), for 1 < i, j < n, 

is AND. Moreover, if n > 2 and no off-diagonal elements of A vanish, then B is 
strictly AND whenever g is strictly AN. 

Proof. 

(1) The matrix A^ = f{\x l — x 3 \ 2 ) is an AND matrix with all diagonal entries 
zero. Hence, by Theorem 2.2.4, we can find n vectors y 1 , . . . , y n G 1Z n such 
that 

But g is a CND1 function, and so the matrix B G 7£ nXn defined by 
B ij = g(\y i -yi\ 2 )=gof(\x i -xi\ 2 ), 
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is also an AND matrix. Thus g o f is a CND1 function. The condition 
that / vanishes only at allows us to deduce that y 1 ^ y- 7 , whenever i ^ j. 
Thus B is strictly AND if g is strictly CND1. 
(2) We observe that A satisfies the hypotheses of Theorem 2.2.4. We may 
therefore write A^ = \y l — y J | 2 , and thus B is AND because g is CND1. 
Now, if A^ 7^ if i ^ j, then the vectors y 1 , ...,y n are distinct, so that B 
is strictly AND if g is strictly CND1. | 

For the next two examples only, we shall need the following concepts. Let 
us call a function g: [0, oo) — > [0, oo) positive definite if, for any positive integers n 
and <i, and for any points x 1 , . . . , x n G 7Z d , the matrix A G 1Z nXn defined by 

A^ = g{\x l - x 3 \ 2 ), forl<z,j<n, 

is non-negative definite. Furthermore, we shall call g strictly positive definite if 
the matrix A is positive definite whenever the points distinct. We 

reiterate that these last two definitions are needed only for examples 2.3.3 and 
2.3.4. 

Example 2.3.3. A Euclidean distance matrix A is AND, indeed strictly so given 
distinct points. This was proved by Schoenberg (1938) and rediscovered by Mic- 
chelli (1986). Schoenberg also proved the stronger result that the matrix 

Aij = \x l — for 1 < i, j < n, 

is strictly AND given distinct points x 1 , . . . , x n G 7Z d , n > 2 and < a < 2. We 
shall derive this fact using Micchelli's methods in Corollary 2.3.7 below, but we 
shall use the result here to illustrate Theorem 2.3.2. We see that, by Theorem 
2.2.4, there exist n vectors y 1 , . . . , y n G 1Z n such that 

A ij = \x i -xi\ a = \y i -yi\ 2 . 

The vectors y 1 , . . . , y n must be distinct whenever the points x 1 , . . . , x n G 1Z d are 
distinct, since A^ ^ whenever i ^ j. 
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Now let g denote any strictly positive definite function. Define B G TZ 



nxn 



by 



B i:i = g(Aij). 



Thus 



g(\x i -xi\ a ) = g(\y i -yi\ 2 ). 



Since we have shown that the vectors y <,..., y n are distinct, the matrix B is 
therefore positive definite. 

For example, the function g(t) = exp(— t) is a strictly positive definite 
function. For an elementary proof of this fact, see Micchelli (1986), p. 15 . Thus 
the matrix whose elements are 



is always (i) non-negative definite, and (ii) positive definite whenever the points 



Example 2.3.4. This will be our first example using a p-norm with p ^ 2. Sup- 
pose we are given distinct points x 1 , . . . , x n G 1Z d . Let us define A G TZ nXn by 



recalling that x\ denotes the k th coordinate of the point x 1 . 
We now remark that A = Ei=i^ (fc) - But 

every 

A (k) 

is a Euclidean dis- 
tance matrix, and so every A^ is AND. Consequently A, being the sum of AND 
matrices, is itself AND. Now A has all diagonal entries zero. Thus, by Theorem 
2.2.4, we can construct n vectors y 1 , . . . , y n G 1Z n such that 



Bij = exp( — \x l — x J | a ), 1 < i, j < n, 



x 1 , . . . , x n are distinct 




x 1 - ar'Hi. 



Furthermore, for k = 1, . . . , d, let A^ G H nXn be given by 





x l — x J ||i = \y l — y 
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As in the preceding example, whenever the points x 1 , . . . ,x n are distinct, so too 
are the vectors y 1 , . . . , y n . 

This does not mean that A is non-singular. Indeed, Dyn, Light and Cheney 
(1991) observe that the 1-norm distance matrix is singular for the distinct points 
{(0,0), (1,0), (1,1), (0,1)}. 

Now let g be any strictly positive definite function. Define B e TZ nXn 

by 

By = g(A tJ ) = g(\\x l - £ J ||i) = g(\y l - y>\ 2 ). 

Thus B is positive definite. 

For example, we see that the matrix By = exp( —\\x l — x J ||i) is positive 
definite whenever the points x 1 , . . . ,x n are distinct. | 

Example 2.3.5. As in the last example, let Ay = \\x l — where n > 2 

and the points x 1 , . . . ,x n are distinct. Now the function /(£) = (1 + t)? is strictly 
CND1 ( Micchelli (1986) ). This is the CND1 function generating the multiquadric 
interpolation matrix. We shall show the matrix B e n nXn defined by 

B ij = f(A ij ) = (l + \\x i -xi\\ 1 )t 

to be strictly AND. 

Firstly, since the points distinct, the previous example shows 

that we may write 

Aij = H^-^'Hi = \y l -y J \ 2 , 

where the vectors y 1 ,...^ 71 are distinct. Thus, since / is strictly CND1, we 
deduce from Definition 2.3.1 that B is a strictly AND matrix. | 

We now return to the main theme of this chapter. Recall that a function 
/ is completely monotonic provided that 

{-l) k f {k \x) > 0, for every k = 0,1,2,... and for < x < oo. 

We now require a theorem of Micchelli (1986), restated in our notation. 
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Theorem 2.3.6. Let f: [0, oo) — > [0, oo) have a completely monotonic derivative. 
Then f is a CND1 function. Further, if f is non-constant, then f is strictly 
CND1. 

Proof. This is Theorem 2.3 of Micchelli (1986). | 

Corollary 2.3.7. The function g(t) = t T is strictly CND1 for every r G (0, 1). 
Proof. The conditions of the previous theorem are satisfied by g. | 

We see now that we may use this choice of g in Theorem 2.3.2, as in the 
following corollary. 

Corollary 2.3.8. For every r G (0,1) and for every positive integer k G [l,d], 
define A™ eTZ nXn by 

Af = \xi-x{\ 2 \ forl<i,j<n. 
Then every A™ i s AND. 

Proof. For each k, the matrix (\x l k — a^.|)" J -_ 1 is a Euclidean distance matrix. 
Using the function g(t) = t T , we now apply Theorem 2.3.2 (2) to deduce that 
A^ = g{\x l - x j \ 2 ) is AND. | 

We shall still use the notation ||.|| p when p G (0, 1), although of course these 
functions are not norms . 

Lemma 2.3.9. For every p G (0, 2), the matrix A G TZ nXn defined by 

A ij = \\ x * ~ xJ \\^ f° r l<i,j <n, 

is AND. If n > 2 and the points x 1 , . . . ,x n are distinct, then we can find distinct 
y 1 , . . . , y n G lZ n such that 

\\ x i- x i\\r p = \yi- y i\\ 

Proof. If we set p = 2r, then we see that r G (0, 1) and A = Yl=i A {k \ where 
the A^ are those matrices defined in Corollary 2.3.8. Hence so that each A^ is 
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AND, and hence so is their sum. Thus, by Theorem 2.2.4, we may write 



A ij = \\x i -xi\\* = \y i -y'\ 2 . 

Furthermore, if n > 2 and the points x 1 , . . . , x n are distinct, then Aij ^ whenever 
i 7^ j, so that the vectors y 1 , . . . , y n are distinct. | 

Corollary 2.3.10. For any p G (0,2) and for any a G (0, 1), define B G TZ nXn 
by 

B ij = (\\x i - x ^r P r. 

Then B is AND. As before, if n > 2 and the points x 1 , . . . , x n are distinct, then 
B is strictly AND. 

Proof. Let A be the matrix of the previous lemma and let g(t) = t T . We now 
apply Theorem 2.3.2 (2) | 

Theorem 2.3.11. For every p G (1,2), the p-norm distance matrix B G 1Z nXn , 
that is: 

Bij = || x 1 - x J || p , for 1 < i, j < n, 

is AND. Moreover, it is strictly AND if n > 2 and the points x 1 ,...,x n are 
distinct, in which case 

(-l) n - 1 detS > 0. 

Proof. If p G (1,2), then a = 1/p g(0,1). Thus we may apply Corollary 
2.3.12. The final inequality follows from the statement of proposition 2.2.3. | 



We may also apply Theorem 2.3.2 to the p— norm distance matrix, for 
p G (1, 2], or indeed to the p th power of the p— norm distance matrix, for p G (0, 2). 
Of course, we do not have a norm for < p < 1, but we define the function in 
the obvious way. We need only note that, in these cases, both classes satisfy 
the conditions of Theorem 2.3.2 (2). We now state this formally for the p— norm 
distance matrix 
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Corollary 2.3.12. Suppose the matrix B is the p—norm distance matrix defined 
in Theorem 2.3.13. Then, if g is a CND1 function, the matrix g(B) defined by 



g(B)ij = g(Bij), for 1 < i, j < n, 

is AND. Further, if n > 2 and the points x 1 ,...,x n are distinct, then g(B) is 
strictly AND whenever g is strictly AN. 

Proof. This is immediate from Theorem 2.3.11 and the statement of Theorem 



2.4. The case p > 2 

We are unable to use the ideas developed in the previous section to understand 
this case. However, numerical experiment suggested the geometry described below, 



TZ m @ TZ n . Given any p > 2, we take the vertices T m of [-m _1 / p , m" 1 ^]™ C TZ m 
and embed this in 1Z m+n . Similarly, we take the vertices Y n of [— n~ 1 ^ p ] n C 
1Z n and embed this too in 1Z m+n . We see that we have constructed two orthogonal 
cubes lying in the p-norm unit sphere. 

Example. If m = 2 and n = 3, then T m = {(±a, ±a, 0, 0, 0)} and T n = 
{(0, 0, ±/3, ±0, ±/3)}, where a = 2" 1 /? and = 3^^. 

Of course, given m and n, we are interested in values of p for which the 
p—norm distance matrix generated by T m U T n is singular. Thus we ask whether 
there exist scalars {\ y }{ ye r m } and {tt z }{zer n }, n °t all zero, such that the function 



vanishes at every interpolation point. In fact, we shall show that there exist scalars 
A and /i, not both zero, for which the function 



2.3.2 (2). 



which proved surprisingly fruitful. We shall view 71 



as two orthogonal slices 



yer m zer„ 
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vanishes at every interpolation point. 
We notice that 
(i) For every y G T m and z G T n , we have \\y — z\\ p = 2 1 / p . 
(ii) The sum J2 y er m \\v ~ v\\p takes the same value for every vertex y G T m , 
and similarly, mutatis mutandis, for T n . 
Thus our interpolation equations reduce to two in number: 

A Yl \\V-V\\p + 2 n+V > = °> 
yer m 

and 

2 m+l/p X + \\z-4p = 0, 

zer n 

where by (ii) above, we see that y and z may be any vertices of F m , T n respectively. 

We now simplify the (1,1) and (2,2) elements of our reduced system by use 
of the following lemma. 

Lemma 2.4.1. Let T denote the vertices of [0, l] k . Then 

K \ ill 
I 



En* = E(;)' 1/p - 

xer i=o 



Proof. Every vertex of V has coordinates taking the values or 1. Thus the 
distinct p- norms occur when exactly / of the coordinates take the value 1, for 
I = 0, . . . , k; each of these occurs with frequency (J . | 

Corollary 2.4.2. 

\\y ~ y Wp = 2 S ( k ) ( k / m ) 1/p ' f° r ever v y e and 

yer m k=o ^ ' 

U-4p = ^( n i )(lM 1/p , for every zeT n . 
zer n i=o ^ ' 

Proof. We simply scale the result of the previous lemma by 2m _1 / p and 2n _1 / p 
respectively. | 
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With this simplification, the matrix of our system becomes 

2Er=o(I)(W 1/p 2".2Vp \ 

2-.2Vp 2£r=o(?)('/») 1/, 7 ' 

We now recall that 

S,(/ p ,1/2)=2-^(^(j7z) 1 / P 
i=o w 

is the Bernstein polynomial approximation of order i to the function f p (t) = t x l p 
at t = 1/2. Our reference for properties for Bernstein polynomial approximation 
will be Davis (1975), sections 6.2 and 6.3. Hence, scaling the determinant of our 
matrix by 2~( m+n \ we obtain the function 

p m , n (p) = 4S m (/ p , l/2)S n (/ p , 1/2) - 2 2 /p. 

We observe that our task reduces to investigation of the zeros of (p m , n - 
We first deal with the case m = n, noting the factorization: 

^ n , n (p) = {2B n (f P , 1/2) + 2 1 ^}{2S n (/ p , 1/2) - 2^}. 

Since / p (t) > 0, for t > we deduce from the monotonicity of the Bernstein 
approximation operator that B n (f p , 1/2) > 0. Thus the zeros of <^ n>n are those of 
the factor 

1> n ( P ) = 2B n (f p ,l/2)-2 1 '*. 

Proposition 2.4.3. ifj n enjoys the following properties. 
(1) ip n (p) — > ip(p), where ip(p) = 2 1_1 / p — 2 1 / p , as n — >■ 00. 

For every p > 1, ip n (p) < ^n+i(p), for every positive integer n. 

(3) For each n, ifj n is strictly increasing for p G [1, 00). 

(4) For every positive integer n, lim^oo ip n (p) = 1 — 2 1_n . 

Proof. 

(1) This is a consequence of the convergence of Bernstein polynomial approxi- 
mation. 
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(2) It suffices to show that B n (f p , 1/2) < B n+ i(f p , 1/2), for p > 1 and n a 
positive integer. We shall use Davis (1975), Theorem 6.3.4: If g is a convex 
function on [0, 1], then B n (g, x) > B n+ i(g, x), for every x G [0, 1]. Further, 
if g is non-linear in each of the intervals [^p, for j = 1, . . . , n, then the 
inequality is strict. Every function f p is concave and non-linear on [0, 1] for 
p > 1, so that this inequality is strict and reversed. 

(3) We recall that 

iPnip) = 2B n (f P , 1/2) - 2 1 /* = J2 (*) - 2VP - 

Now, for p2 > pi > 1, we note that t 1 ^ 2 > t 1 / pi , for i G (0, 1), and also 
that 2 1 / p2 < 2Vpi. Thus (k/n) 1 ^ 2 > (k/n) 1 ^ 1 , for jfe = 1, . . . , n - 1 and so 

Vvfe) > V>n(Pl)- 

(4) We observe that, as p — > oo, 



l = 2(l-2" n )-l = l-2 



l-n 



Corollary 2.4.4. For every integer n > 1, eachifj n has a unique rootp n G (2, oo). 
Further, p n — > 2 strictly monotonically as n — > oo. 

Proof. We first note that ^(2) = 0, and that this is the only root of ip. By 
proposition 2.4.3 (1) and (2), we see that 

lim Vn(2) = ^(2) = and Vn(2) < Vn+i(2) < ^(2) = 0. 

n— >oo 

By proposition 2.4.3 (4), we know that, for n > 1, ip n is positive for all sufficiently 
large p. Since every i(; n is strictly increasing by proposition 2.4.3 (3), we deduce 
that each ip n has a unique root p n G (2, oo) and that ip n (p) < (>)0 for p < (>)p n . 

We now observe that ip n +i(Pn) > i^niPn) = 0, by proposition 2.4.3 (2), 
whence 2 < p n+ i < p n . Thus (p n ) is a monotonic decreasing sequence bounded 
below by 2. Therefore it is convergent with limit in [2,oo). Let p* denote this 
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limit. To prove that p* = 2, it suffices to show that ip(p*) = 0, since 2 is the 
unique root of tp. Now suppose that ip(p*) ^ 0. By continuity, is bounded away 
from zero in some compact neighbourhood N of p* . We now recall the following 
theorem of Dini: If we have a monotonic increasing sequence of continuous real- 
valued functions on a compact metric space with continuous limit function, then 
the convergence is uniform. A proof of this result may be found in many texts, 
for example Hille (1962), p. 78. Thus ty n — > ijj uniformly in N . Hence there is 
an integer no such that ip n is bounded away from zero for every n > uq. But 
p* = limpn and ifj n (Pn) = for each n, so that we have reached a contradiction. 
Therefore ip(p*) = as required. | 

Returning to our original scaled determinant v^n.n, we see that T n U T n 
generates a singular p n -norm distance matrix and p n \ 2 as n — > oo. Furthermore 



using the same method of proof as in proposition 2.4.3 (2). Thus <^ m , n has a unique 
root p m ,n lying hi the interval (p n ,p m ). We have therefore proved the following 
theorem. 

Theorem 2.4.5. For any positive integers m and n, both greater than 1, there is 
a>Pm,n > 2 such that the T m UT n - generated p m ^ n -norm distance matrix is singular. 
Furthermore, if 1 < m < n, then 



and p n \2 as n — > oo. 

Finally, we deal with the "gaps" in the sequence (p n ) as follows. Given 
a positive integer n, we take the configuration T n U r n (0), where T n (9) denotes 
the vertices of the scaled cube [-#n -1 / p , #n -1 / p ] n and 9 > 0. The 2x2 matrix 
deduced from corollary 2.4.2 on page 8 becomes 



P m ,m(p) < <Pm,n(p) < <fn,nip), for 1 < m < n, 



Pm — Pm,m ^ P 



m,n 



Pn,n — Prn 
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Thus, instead of the function (p njTl discussed above, we now consider its analogue: 

<p n , n ,e(p) = 40B 2 M P , 1/2) - (1 + e?) 2/p . 

If p > p n , the unique zero of our original function <^ n>n , we see that <fi n ,n,i(p) = 
<£n,n(p) > 0, because every <^ n>n is strictly increasing, by proposition 2.4.3 (3). 
However, we notice that lim^^o fn,n,e(p) = _ 1, so that ^p n ,n,d(p) < for all 
sufficiently small 9 > 0. Thus there exists a 9* > such that (p n ,n,e* (p) = 0. 
Since this is true for every p > p n , we have strengthened the previous theorem. 
We now state this formally. 

Theorem 2.4.6. For every p > 2, there is a configuration of distinct points 
generating a singular p-norm distance matrix. 

It is interesting to investigate how rapidly the sequence of zeros (p n ) con- 
verges to 2. We shall use Davis (1975), Theorem 6.3.6, which states that, for any 
bounded function / on [0, 1], 

1 

lim n(B n (f,x) — f{x)) = —x(l — x)f"(x), whenever f"(x) exists. 

n— >-oo 2 

Applying this to 

^n(p)=2S n (/ p ,l/2)-2 1 /p, 

we shall derive the following bound. 
Proposition 2.4.7. p n = 2 + 0(n _1 ). 
Proof. We simply note that 

= ^n(Pn) 

= ip(p n ) + 0{n~ l ), by Davis (1975) 6.3.6, 
= V>(2) + (Pn - 2)V'(2) + o(p n - 2) + Oin- 1 ). 
Since ip'(2) ^ 0, we have p n -2 = O^" 1 ). | 
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3.1. Introduction 

In this chapter we use Fourier transform techniques to derive inequalities of the 
form 



where fi is a positive constant and Y^=i Vi = 0- Here we are using the notation of 
the abstract. It can be shown that equation (3.1) implies the bound || A -1 1| 2 < l//z 
(see Chapter 4). Such estimates have been derived in Ball (1989), Narcowich and 
Ward (1990, 1991) and Sun (1990), using a different technique. The author submits 
that the derivation presented here for the Euclidean norm is more perspicuous. 
Further, we relate the generalized Fourier transform to the measure that occurs in 
an important characterization theorem for those functions <p considered here. This 
is useful because tables of generalized Fourier transforms are widely available, thus 
avoiding several of the technical calculations of Narcowich (1990, 1991). Finally, 
we mention some recent work of the author that provides the least upper bound 
on ||A _1 ||2 when the points (xj)j EZ d form a subset of Z d . 

The norm || • || will always be the Euclidean norm in this section. We shall 
denote the inner product of two vectors x and y by xy. 

3.2. The Univariate Case for the Euclidean Norm 

Let n > 2 and let (xj)™ be points in 1Z satisfying the condition \\xj — Xk\\ > 1 for 
j 7^ k. We shall prove that 



whenever Y^=i Vj = 0- 

We shall use the fact that the generalized Fourier transform of (f(x) = \x\ 
is ip(t) = —2/t 2 in the univariate case. A proof of this may be found in Jones 
(1982), Theorem 7.32. 



y T Ay < -\xy T y, 



yeTZ n , 



(3.1) 




j,k=l 



31 



Norm estimates for distance matrices 
Proposition 3.2.1. // Y^=iVj = 0, then 



n poo n 

^2 VjVk \\xj -x k \\ = (2n)~ 1 / (-2/t 2 ) ^ VjVk exp(i(x J - - x k )t) dt 
j.k=i j,k=i 

3.2 

/oo n 
l^y^lH- 2 dt. 
-oo J = 1 



Proof. The two expressions on the righthand side above are equal because of the 
useful identity 

n n 

VjVk exp(i( Xj -x k )t) = \J2vj e ' X3t \ 2 - 

j,k=l j=l 

This identity will be used several times below. We now let 

n 
3 = 1 

The condition Y^=i Vj = implies that g is uniformly bounded. Further, since 
git) = 0(t~ 2 ) for large \t\, we see that g is absolutely integrable. Thus we have 



the equation 



/oo 
g(t) exp(ixt) dt. 
-oo 



A standard result of the theory of generalized Fourier transforms (cf. Jones (1982), 
Theorem 7.14, pages 224ff) provides the expression 

y J yk\\x + x J -x k \\ = (2n)- 1 (-2t- 2 )\J2 yj e ix ^\ 2 exp(ixt)dt, 

3,k=l J -°° 3 = 1 

where we have used the identity stated at the beginning of this proof. We now 
need only set x = in this final equation. | 

Proposition 3.2.2. Let B : 1Z — >■ 1Z be a continuous function such that supp (B) 
is contained in the interval [—1, 1] and < B(t) < t~ 2 . If n > 2, \\xj — Xk\\ > 1 
for j ^ k, and Y!j=i Vj = °> then 

n 

VjVk \\xj -x k \\ < -2B(0)\\y\\ 2 . 

3,k=l 
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Proof. By Proposition 3.2.1 and properties of Fourier transforms, 

n poo n 

y i yk W x 3 ~ x k\\ < ( 27r ) _1 / (-^Bit)) ^2 VjVk exp(i(xj - x k )t)dt 
j,k=i J ~°° j,k=i 

n 

= -2 ^2 VjVkBixj -x k ) 

j,k=l 
= -2B(o)\\y\\ 2 , 

where the first inequality follows from the condition B(t) < t~ 2 . The last line is a 
consequence of supp(S) C [—1, 1]. | 



Corollary 3.2.3. Let 

B(x) 



(l-|x|)/4, i/|x|<l 
0, otherwise. 



Then B satisfies the conditions of Proposition 3.2.2 and B(0) = 1/4. 
Proof. By direct calculation, we find that 

. _ sin 2 (t/2) 1 

m ~ - ¥■ 

It is clear that the other conditions of Proposition 3.2.2 are satisfied. | 
We have therefore shown the following theorem to be true. 

Theorem 3.2.4. Let (xj)™ be points in 1Z such that n > 2 and \\xj — x k \\ > 1 
when j 7^ k. If Y^Jj=i Vj = 0? then 



j,k=l 



We see that a consequence of this result is the non-singularity of the Euclidean 
distance matrix when the points (xj)™ are distinct and n > 2. It is important 
to realise that the homogeneity of the Euclidean norm allows us to replace the 
condition u \\xj — x k \\ > 1 if j ^ k" by u \\xj — x k \\ > e if j ^ k" . We restate 
Theorem 3.2.4 in this form for the convenience of the reader: 
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Theorem 3.2.4b. Choose any e > and let (xj)™ be points in 1Z such that n > 2 
and \\xj — Xk\\ > e when j ^ k. If Y^=i Uj = ®> then 



1 

We shall now show that this bound is optimal. Without loss of generality, we 
return to the case e = 1. We take our points to be the integers 0, 1, . . . , n, so that 
the Euclidean distance matrix, A n say, is given by 



1 



1 




2 
1 



A n — 

\n n — 1 n — 2 

It is straightforward to calculate the inverse of A n : 

f(\-n)/2n 1/2 

1/2 -1 1/2 

1/2 -1 

A- 1 



n \ 
n-1 

J 



l/2n \ 



-1 1/2 
V l/2n 1/2 (l-n)/2nJ 

Proposition 3.2.5. We have the inequality 2 — (it 2 /2n 2 ) < 1 1 -r4.~ 1 1 1 2 < 2. 

Proof. We observe that j | 1 ] 1 2 < ll^-n^li = 2, establishing the upper bound. For 

the lower bound, we focus attention on the (n — l)x(n — 1) symmetric tridiagonal 

minor of A~ x formed by deleting its first and last rows and columns, which we 

shall denote by T n . Thus we have 

(-1 1/2 ^ 
1/2 -1 1/2 
1/2 -1 



V 



-1 1/2 
1/2 -1 / 



Now 



|T n || 2 = max{y T A n 1 y : y T y = 1 and y 1 = y n+1 = 0} 
< max{y T A~ 1 y : y T y = 1} 

= ll^ 1 ||2, 
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so that ||T n || 2 < || A" 1 1| 2 < 2. But the eigenvalues of T n are given by 

Afc = — 1 + cos(/c7r/n), for k = 1, 2, . . . , n — 1. 

Thus ||T n || 2 = 1 — cos(7r — ir/n) > 2 — n 2 /2n 2 , where we have used an elementary 
inequality based on the Taylor series for the cosine function. The proposition is 
proved. | 



3.3. The Multivariate Case for the Euclidean Norm 

We first prove the multivariate versions of Propositions 3.2.1 and 3.2.2, which gen- 
eralize in a very straightforward way. We shall require the fact that the generalized 
Fourier transform of (p(x) = \\x\\ in 1Z d is given by 

0(t) = -c d \\t\\- d -\ 

where 

c d = 2 d 7r( d - 1 )/ 2 r((rf+l)/2). 

This may be found in Jones (1982), Theorem 7.32. We now deal with the analogue 
of Proposition 3.2.1. 

Proposition 3.3.1. // Y^j=i Vj = 0? then 

n „ n 

y J y k \\x J -x k \\=-c d (2n)- d / I^T^e^l 2 !!^-^ 1 ^. (3.3) 

3,k=l 3 = 1 

Proof. We define 

n 

g(t) = -c d \\t\\- d - 1 \Y,y^ Xjt \ 2 - 

i=i 

The condition Yll=i Vj = implies this function is uniformly bounded and the de- 
cay for large argument is sufficient to ensure absolute integrability. The argument 
now follows the proof of Proposition 3.2.1, with obvious minor changes. | 
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Proposition 3.3.2. Let B : 1Z d — > 1Z be a continuous function such that supp[B) 
is contained in the ball {x £ lZ d : < 1}, < B(t) < ||t|r d_1 and B(0) > 0. // 
n > 2, \\xj — Xk\\ > 1 for j ^ k, and Y^j=i Vi = 0? then 

n 

^2 y i yk ~ x k\\ < - c dB(Q)\\y\\ 2 . 

j,k=i 

Proof. The proof of Proposition 3.2.2 clearly generalizes to this case. | 

However, to exhibit a function B satisfying the conditions of Proposition 
3.3.2 is harder than in the univariate case. We modify a construction of Narcowich 
and Ward (1990) and Sun (1990). Let 



B (x) 



1, if||x||<l/2 
0, otherwise. 



Then, using Narcowich and Ward (1990), equation 1.10 or [9], Lemma 3.3.1, we 
find that 

Bo(t) = m\t\\)- d/2 J i (\\t\\m, 

where Jfc denotes the k th -order Bessel function of the first kind. Further, Bq is a 
radially symmetric function since Bq is radially symmetric. We now define 

B = B * S , 

so that, by the convolution theorem, 

B(t) = (B f(t) 

= (2||t||)- d j|(||t||/2), 

2 

and the behaviour of Jo for large argument provides the inequality 

B(t) < lid\\t\\- d -\ 

for some constant fid- Since the conditions of Proposition 3.3.2 are now easy to 
verify when B is scaled by /i^ 1 , we see that we are done . 
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3.4. Fourier Transforms and Bessel Transforms 

Here we relate our technique to the work of Ball (1989) and Narcowich and Ward 
(1990, 1991). 

Definition 3.4.1. A real sequence (yj)j e z d is sa *id to be zero-summing if it is 
finitely supported and J2jez d Vi = 0- 

Definition 3.4.2. A function ip: [0, oo) — > 1Z will be said to be conditionally neg- 
ative definite of order 1 on lZ d , hereafter shortened to CNDl(e£), if it is continuous 
and, for any points (xj)j EZ d in Tt d an d any zero-summing sequence (yj)j e z d i we 
have 

VjVk^Wxj -x k \\) < 0. 

j,kez d 

Such functions were characterized by von Neumann and Schoenberg (1941). For 
every positive integer d, let O d : [0, oo) — > 1Z be defined by 

tod(r) = u d -i~ l / cos(ryu)dy, 

where u may be any unit vector in 1Z d , S^ 1 denotes the unit sphere in 1Z d , and 
ua-i its (d — l)-dimensional Lebesgue measure. Thus fid is essentially the Fourier 
transform of the normalized rotation invariant measure on the unit sphere. 

Theorem 3.4.3. Let ip: [0, oo) — >■ 1Z be a continuous function. A necessary and 
sufficient condition that p be a CNDl(d) function is that it have the form 

/•OO 

ip(r) = <p(0) + (1 - n d (rt)) t~ 2 dP(t), 
Jo 

for every r > 0, where f3: [0, oo) — >■ 1Z is a non- decreasing function such that 
t~ 2 d(3(t) < oo and (3(0) = 0. Furthermore, (3 is uniquely determined by cp. 

Proof. The first part of this result is Theorem 7 of von Neumann and Schoenberg 
(1941), restated in our terminology. The uniqueness of (3 is a consequence of 
Lemma 2 of that paper. | 
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It is a consequence of this theorem that there exist constants A and B such that 
(f(r) < Ar 2 + B. For we have 

/OO />CSO 
(l-n d (rt))t- 2 d(3(t) <2 J t~ 2 d(3(t) < oo, 

using the fact that |0^(r)| < 1 for every r > 0. Further, we see that 
< 1 - n d (p) = 2u d - 1 - 1 / sm 2 (put/2) dt < p 2 /2, 

which provides the bound 

jf 1 (1 - n d (rt)) t~ 2 dP(t) | < r 2 J 1 = 

Thus A = /3(l)/2 and S = <p(0) + 2 t~ 2 dfi(t) suffice. Therefore the function 
: x G lZ d } is a tempered distribution in the sense of Schwartz (1966) and 
possesses a generalized Fourier transform {y?(||£||) : £ G 1Z d }. There is a rather 
simple relation between the generalized Fourier transform and the nondecreasing 
function of Theorem 3.4.3 for a certain class of functions. This is our next topic. 

Definition 3.4.4. A function <p: [0, oo) — > TZ will be termed admissible if it is a 
continuous function of algebraic growth which satisfies the following conditions: 

1. is a continuous function on lZ d \ {0}. 

2. The limit\im m ^ ||£|| d+ V(ll£ll) exists. 

3. The integral /{||£||>i} l^dlClDI ^£ exists. 

It is straightforward to prove the analogue of Propositions 3.2.1 and 3.3.1 
for an admissible function. 

Proposition 3.4.5. Let <p: [0, oo) — >■ 1Z be an admissible function and let (yj)j e z d 
be a zero-summing sequence. Then for any choice of points (xj)j EZ d in lZ d we have 
the identity 

E yjvMWxj - x k \\) = W / J E ^exp^ol^dlell)^- (3.4) 

j,kez d nd jez d 
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Proof. Let g: 1Z d — >■ 1Z be the function denned by 

l 2 

Then g is an absolutely integrable function on TZ d , because of the conditions on 
ip and because {%ij)j^z d is a zero-summing sequence. Thus g is the generalized 
transform of fc j/jj/^y ( || • +£j — ^fcll)) an d by standard properties of generalized 
Fourier transforms we deduce that 

^2yjVk<p(\\x + Xj -x k \\) = (2n)- d [ I Vjexpfajg) <p(\\£\\) exp(ix£) d£. 

The proof is completed by setting x = 0. | 
Proposition 3.4.6. Let (p: [0, oo) 71 be an admissible CNDl(d) function. Then 

df3(t) = -(27r)-%_ 1 ^(||^||)t d+1 dt, 
where u may be any unit vector in lZ d . 

Proof. Let \i and v be different integers and let (yj)j e z d be a sequence with only 
two nonzero elements, namely y^ = —y v = 2 -1 / 2 . Choose any point ( G 7Z d and 
set x M = 0, x v = C, so that equation (3.4) provides the expression 

^(o) - pencil) = (27r)- d / (i - cos(co) ^(iieii) dC 

Jn d 

Employing spherical polar coordinates, this integral takes the form 

poo 

m - ^(iidl) = 27r)- d Wd _ 1 / (i - n d (f ||CID) ^UHI)^ -1 

where -u may be any unit vector in 7Z d . Setting r = ||£||, we have 

/•OO 

<p(r) = <p{0) + (1 - Q d (rt)) 7 (t)t" 2 d*, 
Jo 

where 7(f) = — (2n)~ d ujd-i0(\\tu\\)t d+1 . Now Theorem 4.2.6 of the following chap- 
ter implies that <p is a nonpositive function. Thus there exists a nondecreasing 
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function (3: [0, oo) — > 1Z such that ^(t) dt = d/3(t), and t~ 2 d(3{i) is finite and 
/3(0) = 0. But the uniqueness of the representation of Theorem 3.4.3 implies that 
(3 = j3, that is 

dp(t) = -(2n)- d u d - 1 £(\\tu\\)t d+1 dt, 
and the proof is complete. | 

This proposition is useful if we want to calculate (3 for a particular function 
(p, since tables of generalized Fourier transforms are readily available. 

Example 3.4.7. Letcp(r) = (r 2 + l) 1/2 . This is a non-negative CNDl(d) function 
for all d (see Micchelli (1986)). When d = 3, the generalized Fourier transform is 
ip(r) = — Aitr~ 2 K 2 (r) . Here K 2 is a modified Bessel function which is positive and 
smooth in 7Z + , has a pole at the origin, and decays exponentially (See Abramowitz 
and Stegun (1970)). Consequently (p is a non-negative admissible function. Ap- 
plying Theorem 3.4-7 gives the equation 

dfi(r) = (2n)- 3 (4n)r 4 (4nr- 2 K 2 (r)) 
= (2r 2 /Ti)K 2 {r)dr, 
agreeing with Narcowich and Ward (1991), equation 3.12. 

3.5. The Least Upper Bound for Subsets of the Integer Grid 

In the next chapter we use extensions of the technique provided here to derive the 
the following result. 

Theorem 3.5.1. Let ip: [0, oo) — > TZ be an admissible function that is not identi- 
cally zero, letip(Q) > 0, and let ip be CND 1(d) for every positive integer d. Further, 
let (xj)j eZ d be any elements of Z d and let A = (<p>(\\xj —Xk\\))j,keAf> where Af can 
be any finite subset of Z d . Then we have the inequality 

P -1 ||<(£ Mhe + 2nk\\)\y\ 

kez d 

where e = [1,...,1] T G lZ d and (p is the generalized Fourier transform of ip. 
Moreover, this is the least upper bound valid for all finite subsets of Z d . 
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Proof. See Section 4.4 of the thesis. 
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4 : Norm estimates for Toeplitz distance matrices I 
4.1. Introduction 

The multivariate interpolation problem is as follows: given points (x_,-)" =1 in TZ d 
and real numbers (/j)j = i, construct a function s: lZ d — > 1Z such that s(xk) = fk, 
for k = l,...,n. The radial basis function approach is to choose a univariate 
function <p: [0, oo) — > 1Z, a norm || . || on lZ d , and to let s take the form 

n 

s ( x ) = ^2vM\\ x - x j\\)- 

i=i 

The norm || . || will be the Euclidean norm throughout this chapter. Thus the 
radial basis function interpolation problem has a unique solution for any given 
scalars (/j)7=i ^ an< ^ oru y ^ the matrix (<p(||x_j — #fc|| ))?*,=! is invertible. Such 
a matrix will, as before, be called a distance matrix. These functions provide a 
useful and flexible form for multivariate approximation, but their approximation 
power as a space of functions is not addressed here. 

A powerful and elegant theory was developed by I. J. Schoenberg and oth- 
ers some fifty years ago which may be used to analyse the singularity of distance 
matrices. Indeed, in Schoenberg (1938) it was shown that the Euclidean dis- 
tance matrix, which is the case cp(r) = r, is invertible if n > 2 and the points 
(^j)j=i are distinct. Further, extensions of this work by Micchelli (1986) proved 
that the distance matrix is invertible for several classes of functions, including the 
Hardy multiquadric, the only restrictions on the points {xj)^ =1 being that they 
are distinct and that n > 2. Thus the singularity of the distance matrix has been 
successfully investigated for many useful radial basis functions. In this chapter, we 
bound the eigenvalue of smallest modulus for certain distance matrices. Specif- 
ically, we provide the greatest lower bound on the moduli of the eigenvalues in 
the case when the points {xj)^ =1 form a subset of the integers Z d , our method 
of analysis applying to a wide class of functions which includes the multiquadric. 
More precisely, let iV be any finite subset of the integers Z d and let A^ in be the 
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smallest eigenvalue in modulus of the distance matrix ((p(\\j — k\\))j : keN- Then 
the results of Sections 3 and 4 provide the inequality 

|A^ in | > CV, (4.1.1) 

where is a positive constant for which an elegant formula is derived. We also 
provide a constructive proof that C v cannot be replaced by any larger number, 
and it is for this reason that we shall describe inequality (4.1.1) as an optimal 
lower bound. Similarly, we shall say that an upper bound is optimal if none of the 
constants appearing in the inequality can be replaced by smaller numbers. 

It is crucial to our analysis that the distance matrix (<p(\\j — k\\))j t keN may 
be embedded in the bi-infmite matrix (<p(\\j ' — k\\))j,kez d - Such a bi-infinite matrix 
is called a Toeplitz matrix if d = 1. We shall use this name for all values of d, 
since we use the multivariate form of the Fourier analysis of Toeplitz forms (see 
Grenander and Szego (1984)). 

Of course, inequality (4.1.1) also provides an upper bound on the norm 
of the inverse of the distance matrices generated by finite subsets of the integers 
Z d . This is not the first paper to address the problem of bounding the norms 
of inverses of distance matrices and we acknowledge the papers of Ball (1989) 
and Narcowich and Ward [1990, 1991], which first interested the author in such 
estimates. Their results are not limited to the case when the data points are a 
subset of the integers. Instead, they apply when the points satisfy the condition 
ll^j — %k\\ > e f° r 3 7^ k, where e is a positive constant, and they provide lower 
bounds on the smallest modulus of an eigenvalue for several functions <p, including 
the multiquadric. We will find that these bounds are not optimal, except in the 
special case of the Euclidean norm in the univariate case. Further, our bounds 
apply to all the conditionally negative definite functions of order 1. The definition 
of this class of functions may be found in Section 4.3. 

As in the previous section, we make extensive use of the theory of gen- 
eralized Fourier transforms, for which our principal reference will still be Jones 
(1982). These transforms are precisely the Fourier transforms of tempered distri- 
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butions constructed in Schwartz (1966). First, however, Section 2 presents several 
theorems which require only the classical theory of the Fourier transform. These 
results will be necessary in Section 4.3. 

4.2. Toeplitz forms and Theta functions 

We require several properties of the Fejer kernel, which is defined as follows. For 
each positive integer n, the n th univariate Fejer kernel is the positive trigonometric 
polynomial 

n 

K n (t) = (! ~ \k\/n)exp(ikt) 



k-- 



(4.2.1) 



sin nt/2 
n sin 2 t/2 

Further, the n th multivariate Fejer kernel is defined by the product 

K n (t 1 ,...,t d ) = K n (t 1 )K n (t 2 )---K n (t d ), ten d . (4.2.2) 

Lemma 4.2.1. The univariate kernel enjoys the following property: for any con- 
tinuous 2rr -periodic function f:TZ-^7l and for all x G 1Z we have 

r-2-K 

lim (27T)" 1 / K n {t - x)f(t) dt = f{x). 

Moreover, we have the equations 

(2/t)- 1 / K n (t)dt = 1 (4.2.3) 
Jo 

and 

n-i 2 

K n (t)= n- 1/2 J2exp(ikt) . (4.2.4) 

fc=0 

Proof. Most text-books on harmonic analysis contain the first property and (4.2.3). 
For example, see pages 89ff, volume I, Zygmund (1979). It is elementary to deduce 
(4.2.4) from (4.2.1). | 

Lemma 4.2.2. For every continuous [0, 27r] d -periodic function f:7Z d — > 1Z, the 
multivariate Fejer kernel gives the convergence property 

lim (27r)" d / K n (t - x)f(t) dt = fix) 

n ^°° J[0,27T]d 
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for every x G 1Z d . Further, K n is the square of the modulus of a trigonometric 
polynomial with real coefficients and 



K n (t)dt = 1. 



Proof. The first property is Theorem 1.20 of chapter 17 of Zygmund (1979). The 
last part of the lemma is an immediate consequence of (4.2.3), (4.2.4) and the 
definition of the multivariate Fejer kernel. | 

All sequences will be real sequences here. Further, we shall say that a 
sequence (o,j) Z d := {aj}j EZ d is finitely supported if it contains only finitely many 
nonzero terms. The scalar product of two vectors x and y in lZ d will be denoted 
by xy. 

Proposition 4.2.3. Let f: lZ d — >■ 1Z be an absolutely integrable continuous func- 
tion whose Fourier transform f is also absolutely integrable. Then for any finitely 
supported sequence (aj) Z d, and for any choice of points (xj) Z d in lZ d , we have the 
identity 

X a 3 a k f{xj - x k ) = (27r) _d ! J ^ aj exp(za^) /(£) d£,. 
j,kez d nd jez d 

Proof. The function ■ k ajakf(x + Xj — Xk) '■ x G lZ d } is absolutely integrable. 
Its Fourier transform is given by 

X ajOkfi' + Xj-Xk) (0= X a J a fc exp(z(x j - x k )0f(0 
j,kez d j,kez d 

= | a.exp^ofim £ 

jez d 

and is therefore absolutely integrable. Therefore the Fourier inversion theorem 
states that 

X ajOkfix + Xj - x k ) = (27r)" d I | X aj exp(ixj£) /(£) exp(z^) d£. 
j,kez d nd jez d 

Setting x = produces the stated equation. | 
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In this dissertation a key role will be played by the symbol function 

*(0= E /(£ + 2t*0, (4.2.5) 

fce.z d 

If / e L^iTV 1 ), then a is an absolutely integrable function on [0, 27r] d and its 
defining series is absolutely convergent almost everywhere. These facts are conse- 
quences of the relations 



oo 



> / 1/(01 dt= E / i/(e+27rfc)|de= / e \m+^k)\dt, 

JK d keZd J[0,27r]d J[0,27r]d keZd 

the exchange of integration and summation being a consequence of Fubini's theo- 
rem. If the points {xj)z d are integers, then we readily deduce the following bounds 
on the quadratic form. 

Proposition 4.2.4. Let f satisfy the conditions of Proposition ^.2.3 and let 

(a,j) Z d be a finitely supported sequence. Then we have the identity 

E a J a k f(j-k) = (2n)- d [ I £ a, exp(u£) d£. (4.2.6) 

Further, letting m = inf{a(0 : £ G [0, 2n] d } and M = sup{a(0 : f G [0, 2yr] d }, we 
/mi>e £/ie bounds 

m E a ? ^ E a jak f(j -k)<M 4 

je.z d j,kez d jez d 

Proof. Proposition 4.2.3 implies the equation 

E a J a kfU-k) = E ( 2yr )" d / I E a J exp(zjO| 2 /(^ + 27rfc)^ 



jez d 

the exchange of integration and summation being justified by Fubini's theorem. 
For the upper bound, the Parseval theorem yields the expressions 

E a 3 a kfU ~k) = {2-K)~ d [ I % exp(zj'0 <r(f) ^ 

< M E 4 
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The lower bound follows similarly and the proof is complete. 



The inequalities of the last proposition enjoy the following optimality prop- 
erty. 

Proposition 4.2.5. Let f satisfy the conditions of Proposition 4-2.3 and suppose 
that the symbol function is continuous. Then the inequalities of Proposition 4-2-4 
are optimal lower and upper bounds. 

Proof. Let £m G [0,2-7r] d be a point such that ct(£m) = M, which exists by con- 
tinuity of the symbol function. We shall construct finitely supported sequences 
{(a>j^)jez d '■ n = 1,2, . . .} such that ^j eZ d{a^) 2 = 1, for all n, and 



lim "j % fti~k) = M. (4.2.7) 

j,kez d 



We recall from Lemma 4.2.2 that the multivariate Fejer kernel is the square 
of the modulus of a trigonometric polynomial with real coefficients. Therefore 
there exists a finitely supported sequence {a!j l ')z d satisfying the relation 

\J2 4 n) ex P (^e)| 2 =^n(e-eM), ien d . (4.2.8) 

jez d 

Further, the Parseval theorem and Lemma 4.2.2 provide the equations 

(af) 2 = (2n)~ d [ K n (£ - £ M ) d£ = 1 

jeZ d J[0,27r] d 

and 

lim (2yr)- d [ K n {i - £mMO di = a(t M ) = M. 



'[0,2n] d 

It follows from (4.2.6) and (4.2.8) that the limit (4.2.7) holds. The lower bound 
of Proposition 4.2.4 is dealt with in the same fashion. | 

The set of functions satisfying the conditions of Proposition 4.2.5 is non- 
void. For example, suppose that we have /(£) = C(||^|| _d ~ 5 ), for large ||£||, where 
5 is a positive constant. Then the series defining the symbol function a converges 
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uniformly, by the Weierstrass M-test, and a is continuous, being a uniformly con- 
vergent sum of continuous functions. These remarks apply when / is a Gaussian, 
which is the subject of the rest of this section. We shall see that the analysis of 
the Gaussian provides the key to many of our results. 

Proposition 4.2.6. Let X be a positive constant and let f(x) = exp ( — A 1 1 x 1 1 2 ) , for 

x G lZ d . Then f satisfies the conditions of Proposition 4-2.5. 

Proof The Fourier transform of / is the function /(£) = (n/X) d / 2 exp(-||£|| 2 /4A), 
which is a standard calculation of the classical theory of the Fourier transform. It 
is clear that / satisfies the conditions of Proposition 4.2.3, and that the symbol 
function is the expression 

*(0 = (^M) d/2 E e M~U + 2^|| 2 /4A), £ G U d . (4.2.9) 

kez d 

Finally, the decay of the Gaussian ensures that a is continuous, being a uniformly 
convergent sum of continuous functions. | 

This result is of little use unless we know the minimum and maximum 
values of the symbol function for the Gaussian. Therefore we show next that 
explicit expressions for these numbers may be calculated from properties of Theta 
functions. Lemmata 4.2.7 and 4.2.8 address the cases when d = 1 and d > 1 
respectively. 

Lemma 4.2.7. Let X be a positive constant and let E\ : 1Z — > 1Z be the 2tt -periodic 
function 

oo 

E 1 (t)= E exp (-A(t + 2/ctt) 2 ) . 

fc= — oo 

Then E^O) > E^t) > Ei(n) for all t G K. 

Proof. An application of the Poisson summation formula provides the relation 

oo 

E x (t) = (4ttA)- 1 /2 e~ k2 / 4X e M 

k= — oo 

= (4ttA)- 1 / 2 ^l + 2Ee" fc2/4A cos(H)j . 
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This is a Theta function. Indeed, using the notation of Whittaker and Watson 
(1927), Section 21.11, it is a Theta function of Jacobi type 

oo 

•d${z,q) = 1 + 2^^ cos (2 £^), 
k=i 

where q G C and \q\ < 1. Choosing q = e _1 / 4A we obtain the relation 

E 1 {t) = {A-K\)- l 'Hz{t/2,q). 

The useful product formula for $3: 

00 

3 (z, q) = G Y[ (1 + 2q 2k ~ 1 cos 2z + q 4k ~ 2 ), 
k=i 

where G = UkLii 1 -Q 2k )i is g iven in Whittaker and Watson (1927), Sections 21.3 
and 21.42. Thus 

00 

E 1 (t) = (47rA)- 1 / 2 G JJ(1 + 2q 2k ~ 1 cos t + q" k ~ 2 ), t G K. 
k=i 

Now each term of the infinite product is a decreasing function on the interval 
[0, 7r], which implies that E\ is a decreasing function on [0, n}. Since E\ is an even 
27r-periodic function, we deduce that E\ attains its global minimum at t = tv and 
its maximum at t = 0. ■ 

Lemma 4.2.8. Let X be a positive constant and let E d :1Z d — > 1Z d be the [0,27r] d - 
periodic function given by 

E d (x)= exp(-A||x + 2kn\\ 2 ). 

kez d 

Then E d (0) > E d (x) > E d (ne), where e = [1, 1, . . . , 1] T 

Proof. The key observation is the equation 

d 



E d (x) = Y[E 1 (x k ). 



k=i 
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Thus E d (0) = nLi^i(O) > nti^i(^) = E d (x) > ULiEi(n) = E d (7re), 
using the previous lemma. | 

These lemmata imply that in the Gaussian case the maximum and minimum 
values of the symbol function occur at £ = and £ = 7re respectively, where 
e = [1, . . ., 1] T . Therefore we deduce from formula (4.2.9) that the constants of 
Proposition 4.2.4 are the expressions 

m=(n/X) d/2 exp(-||7re + 27rA;|| 2 /4A) and 

keZd (4.2.10) 
M=(n/X) d / 2 exp(-||7rfc|| 2 /A). 

kez d 



4.3. Conditionally negative definite functions of order 1 

In this section we derive the optimal lower bound on the eigenvalue moduli of the 
distance matrices generated by the integers for a class of functions including the 
Hardy multiquadric. 

Definition 4.3.1. A real sequence {y 3 )z d is sa -id to be zero-summing if it is finitely 
supported and J2jez d Vi = 0- 

Let <p: [0, oo) — > 1Z be a continuous function of algebraic growth. Thus it is 
meaningful to speak of the generalized Fourier transform of the radially symmetric 
function {<p(||a;||) : x G 1Z d }. We denote this transform by {<£(||£||) : £ G TZ d } : 
so emphasizing that it is a radially symmetric distribution, but we note that (p 
depends on d. We shall restrict attention to the collection of functions described 
below. 

Definition 4.3.2. A function (p: [0, oo) — > 1Z will be termed admissible if it is a 
continuous function of algebraic growth which satisfies the following conditions: 

1. <p is a continuous function on lZ d \ {0}. 

2. The limit lim^| Ko exists. 

3. The integral / { ^||> 1} |^(||C||)| d£ exists. 
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It is straightforward to prove the analogue of Proposition 4.2.3 for an ad- 
missible function. 

Proposition 4.3.3. Let <p: [0, oo) — >■ 1Z be an admissible function and let {i)j)z d 
be a zero-summing sequence. Then for any choice of points (xj) z a in lZ d we have 
the identity 

J2 yjyMWxj - x k \\) = (2n)~ d f J J2 Vj<xp(™&\ 2 <P(U\\)<%- (4-3.1) 

j,kez d nd jez d 

Proof. Let g: lZ d — > 1Z be the function defined by 



2 



9® = | Yl 2/i ex P(^jO 0(U\\)- 

jez d 

Then g is an absolutely integrable function on 1Z d , because of the conditions on ip 
and because (yj)z d is a zero-summing sequence. Thus g is the generalized trans- 
form of J2j k yjUk^dl'+Xj—Xk ||), and by standard properties of generalized Fourier 
transforms we deduce that 

^2yjy k <p(\\x + Xj -x k \\) = (27r) _d / I ^ y i exp(ixjf) <^(||CII) exp(iar^) d^. 
j,k Jnd jez d 

The proof is completed by setting x = 0. | 

We come now to the subject that is given in the title of this section. 

Definition 4.3.4. Let (p: [0, oo) — > 1Z be a continuous function. We shall say that 
p is conditionally negative definite of order 1 on every lZ d , hereafter shortened to 
CND1, if we have the inequality 

VjywiWxj ~ x k \\) < 0, 

j,kez d 

for every positive integer d, for every zero-summing sequence (yj)z d an d for any 
choice of points (xj) z <t in 1Z d . 



Such functions were completely characterized by I. J. Schoenberg (1938). 
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Theorem 4.3.5. A continuous function <p: [0, oo) — > 1Z is CND1 if and only if 
there exists a nondecreasing function a: [0, oo) — > 1Z such that 

/•OO 

<p(r) = <p(0)+ [1 -exp(-tr 2 )]t _1 da(t), for r > 0, 
Jo 

and the integral t~ x da(t) exists. 

Proof This is Theorem 6 of Schoenberg (1938). | 
Thus da is a positive Borel measure such that 

r- 1 />oo 

/ da(t) < oo and / t _1 da(t) < oo. 

Further, it is a consequence of this theorem that there exist constants A and B 
such that (p(r) < Ar 2 + B, where A and B are constants. In order to prove this 
assertion we note the elementary inequalities 

/OO POO 
[1 - exp(-tr 2 )]t _1 da(t) < J t' 1 da(t) < oo, 



and 



/ [1 - exp(-tr 2 )]t _1 da(t) < r 2 / da(t). 
Jo Jo 



Thus A = r 2 (a(l) - a(0)) and B = <p(0) + t' 1 da(t) suffice. Therefore we may 
regard a CND1 function as a tempered distribution and it possesses a generalized 
Fourier transform. The following relation between the transform and the integral 
representation of Theorem 4.3.5 will be essential to our needs. 

Theorem 4.3.6. Let ip\ [0, oo) 71 be an admissible CND1 function. For £ G 
lZ d \ {0}, we have the formula 

poo 

m\\) = - / eM-m\ 2 /m*/t) d/2 t- 1 m*)- (4-3.2) 

Before embarking on the proof of this theorem, we require some ground- 
work. We shall say that a function /: 1Z d \ {0} — > 1Z is symmetric if f(—x) = f(x), 
for every x G lZ d \ {0}. 
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Lemma 4.3.7. Let a: [0, oo) — >■ 1Z be a nondecreasing function such that the 
integral t~ l da(t) exists. Then the function 

pea 

HO = - eM-Ufmin/t^t-Uait), £ G lZ d \ {0}, (4.3.3) 
Jo 

is a symmetric smooth function, that is every derivative exists. 
Proof. For every nonzero £, the limit 

limexp(-||^|| 2 /4t)(7rA) d /V 1 = 

implies that the integrand of expression (4.3.3) is a continuous function on [0, oo). 
Therefore it follows from the inequality 

/OO />CSO 
exp(-||^|| 2 /4t)(7r/t) d/2 t- 1 da(t) < n d/2 J t~ l da(t) < oo 

that the integral is well-defined. Further, a similar argument for nonzero £ shows 
that every derivative of the integrand with respect to £ is also absolutely inte- 
grable for t e [0, oo), which implies that every derivative of ip exists. The proof is 
complete, the symmetry of tp being obvious. | 

Lemma 4.3.8. Let f: lZ d — >■ 1Z be a symmetric absolutely integrable function such 
that 



/ a i ex P(^ x j^) f(t) dt = 

Jn d 1 ~~t r , 



jez d 

for every finitely supported sequence {a,j)z d and for any choice of points (xj) Z d. 
Then f must vanish almost everywhere. 

Proof. The given conditions on / imply that the Fourier transform / is a symmetric 
function that satisfies the equation 

^2 aja k f(xj - x k ) = 0, 
j,kez d 

for every finitely supported sequence {aj)z d an d for all points [xj) Z d in lZ d . Let 
a and (3 be different integers and let a a and ap be the only nonzero elements of 
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(a,j) Z d. We now choose any point £ G 1Z d \ {0} and set x a = 0, xp = £, which 
provides the equation 

Therefore /(0) = /(£) = 0, and since £ was arbitrary, / can only be the zero 
function. Consequently / must vanish almost everywhere. | 

Corollary 4.3.9. Let g: lZ d \ {0} — >■ 1Z be a symmetric continuous function such 
that 

[ \j2vj exp(ix^) 2 \g(0\ <% < oo (4.3.4) 

and 

[ exp(^e) ^(0 rf£ = 0, (4.3.5) 

for every zero-summing sequence (yj) Z d and for any choice of points {xj) Z d. Then 
0(0 = for every CeTZ d \{0}. 

Proof. For any integer k G {1, . . . , d} and for any positive real number A, let h be 
the symmetric function 

/i(0=^)sin 2 A£ fc , £e1l d \{0}. 

The relation 

1 1 2 

HO = 9(0 g exp(iA£ fc ) - - exp(-zA£ fc ) 

and condition (4.3.4) imply that h is absolutely integrable. 

Let (a,j) Z d be any real finitely supported sequence and let (bj) Z d be any 
sequence of points in lZ d . We define a real sequence (yj) Z d and points (xj) Z d in 
7£ d by the equation 

^2 yj-exp(ixjO = sinA£ fc a j exp(z6 i £). 
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Thus (yj)z d is a sequence of finite support. Further, setting £ = 0, we deduce that 
J2j eZ dUj = 0, so (yj)z d is a zero-summing sequence. By condition (4.3.5), we 
have 




Therefore we can apply Lemma 4.3.8 to h, finding that it vanishes almost every- 
where. Hence the continuity of g for nonzero argument implies that g(£) sin 2 A£fc = 
for £ 7^ 0. But for every nonzero £ there exist k G {1, . . . , d} and A > such that 
sinA£/c 7^ 0. Consequently g vanishes on 1Z d \ {0}. | 

We now complete the proof of Theorem 4.3.6. 

Proof of Theorem 4-3.6. Let {yj)z d be a zero-summing sequence and let {xj) Z d 
be any set of points in 1Z d . Then Theorem 4.3.5 provides the expression 

^2 VjVkViWxj - x k \\) = - / yjykGxp(-t\\xj - XkW 2 )^- 1 da(t), 

j,kez d Jo j,kez d 

this integral being well-defined because of the condition Yljez d Vj = Therefore, 
using Proposition 4.2.3 with /(•) = exp(— 1\\ ■ || 2 ) in order to restate the Gaussian 
quadratic form in the integrand, we find the equation 



J2 VjVkViWxj -Sfcll) 
j,kez d 




where we have used Fubini's theorem to exchange the order of integration and 
where if) is the function defined in (4.3.3). By comparing this equation with the 
assertion of Proposition 4.3.3, we see that the difference g(£) = <£(||£||) — ^(0 
satisfies the conditions of Corollary 4.3.9. Hence <£(||£||) = ip(£) for all £ e lZ d \ {0}. 
The proof is complete. | 
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Remark. An immediate consequence of this theorem is that the generalized 
Fourier transform of an admissible CND1 function cannot change sign. 

The appearance of the Gaussian quadratic form in the proof of Theorem 
4.3.6 enables us to use the bounds of Lemma 4.2.8, which gives the following result. 

Theorem 4.3.10. Let <p: [0, oo) — >■ 1Z be an admissible CND1 function and let 
(Vj)z d be a zero- summing sequence. Then we have the inequality 



j,kez d jez d 



where e = [!,...,!] 



Proof. Applying (4.3.1) and dissecting 1Z d into integer translates of [0, 2ix] d , we 
obtain the equations 

j,kez d Jnd jez d 

where the interchange of summation and integration is justified by Fubini's the- 
orem, and where we have used the fact that (p does not change sign. Here the 
symbol function has the usual form (4.2.5). Further, using (4.3.2), we again apply 
Fubini's theorem to deduce the formula 

koi= E i^iie +2^11)1 

kez d 

= r(j2 e M-U + ^kf/4t)) (n/tf/^dait). 
It follows from Lemma 4.2.8 that we have the bound 

»oo 



a(OI > / (E exp(-||7re + 27r£;|| 2 /4t)) (tt /t^H^dait) 



kez d 
|cr(7re) | 



(4.3.7) 



The required inequality is now a consequence of (4.3.6) and the Parseval relation 
(2n)~ d [ I £ Vj exp(tj0 V = £ y). 



jez d 
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When the symbol function is continuous on 1Z d \ 2nZ d , we can show that 
the previous inequality is optimal using a modification of the proof of Proposition 

,<->: 

sequences such that lim n ^ 00 J^jezdiVj™') 2 = 1 an d 



4.2.5. Specifically, we construct a set {(yj ') Z d : n = 1,2,...} of zero-summing 

,(»)\2 



I™ I £ ^ n V(lb--fc||) = k( 7 re)|, 



which implies that we cannot replace |cr(7re)| by any larger number in Theorem 
4.3.10. 

Corollary 4.3.11. Let ip: [0, oo) — > 1Z satisfy the conditions of Theorem 4-3.10 
and let the symbol function be continuous in the set lZ d \ 2ivZ d . Then the bound 
of Theorem 4-3.10 is optimal. 

Proof. Let m be an integer such that 4m > d + 1 and let S m be the trigonometric 
polynomial 

d 

^(o = [rf- 1 E sin2 fe/ 2 )] 2m ' t end - 

Recalling from Lemma 4.2.2 that the multivariate Fejer kernel is the square of the 
modulus of a trigonometric polynomial with real coefficients, we choose a finitely 
supported sequence (yY^z- 1 satisfying the equations 

| Y, Vj n) exp^)f = K n{t, ~ Ke)S m (0, * G K d . (4.3.8) 

jez d 

Further, setting £ = we see that (Vj )z d is a zero-summing sequence. Applying 
(4.3.6), we find the relation 

£ 2/j n) ^ n V(lb' - k W) = W d I K n(t ~ ™)S m (t) k(OI (4-3-9) 

Moreover, because the second condition of Definition 4.3.2 implies that <!? m |<7| is a 
continuous function, Lemma 4.2.2 provides the equations 



lim (27v)- d / K n (i - ne)S m (0 \a(0\ d£ = S m (ne) \a(ire)\ = 

'[0,27r] d 
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It follows from (4.3.9) that we have the limit 



l T\ E yf ) y { k ) v{\\3-k\\) = \a^e) 

n— >-oo| J 
j,keZ d 



Finally, since S m is a continuous function, another application of Lemma 
4.2.2 yields the equation 

hm (2n)- d ! K n (t - 7re)S m (£) <% = S m (ive) = 1. 

n ^°° J[0,2n]* 

By substituting expression (4.3.8) into the left hand side and employing the Par- 
seval relation 



(2-)" d / I E yj n) 2 * = E (^ n) ) 

we find the relation lim n ^ QO Zlj e ^d(yj n) ) 2 = 1. 



4.4. Applications 

This section relates the optimal inequality given in Theorem 4.3.10 to the spec- 
trum of the distance matrix, using an approach due to Ball (1989). We apply the 
following theorem. 

Theorem 4.4.1. Let A e TZ nXn be a symmetric matrix with eigenvalues Ai > 
• • • > A n . Let E be any subspace oflZ n of dimension m . Then we have the inequality 

max{a; T AE : x T x = 1, x _L E} > X m +i- 

Proof. This is the Courant-Fischer minimax theorem. See Wilkinson (1965), pages 
99ff. ■ 

For any finite subset N of Z d , let An be the distance matrix (v?(||j — 
k\\))j,keN- Further, let the eigenvalues of A N be Ai > • • • > \\n\, where \N\ is the 
cardinality of N, and let A^ in be the smallest eigenvalue in modulus. 
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Proposition 4.4.2. Let (p: [0, oo) — >■ 1Z be a CND1 function that is not identically 
zero. Let <p(Q) > and let fi be a positive constant such that 

vjvM\\j-k\\)<-n yl ( 4A1 ) 

j,kez d jez d 
for every zero- summing sequence (yj)z d - Then for every finite subset N of Z d we 
have the bound 

l^minl > A*- 
Proof. Equation (4.4.1) implies that 

y T A N y < -\iy T y, 

for every vector (yj)j eN such that YljeNVj = ®- Thus Theorem 4.4.1 implies that 
the eigenvalues of An satisfy —ji > A2 > • • • > \\n\, where the subspace E of 
that theorem is simply the span of the vector [1, 1, ... , 1] T G 1Z N . In particular, 
> A 2 > • • • > A|jv|. This observation and the condition ip(0) > provide the 
expressions 

\N\ \N\ 

< trace A N = Ai + ^ Xj = Ai — ^ \ Xj\. 

3=2 3=2 

Hence we have the relations A^ in = A2 < —fi. The proof is complete. | 

We now turn to the case of the multiquadric <p c (r) = (r 2 + c 2 ) 1 / 2 , in or- 
der to furnish a practical example of the above theory. This is a non-negative 
CND1 function (see Micchelli (1986)) and its generalized Fourier transform is the 
expression 

^c(iieii) = -^-\^c/\m {d+1)/2 K {d+1)/ M\m, 

for nonzero £, which may be found in Jones (1982). Here {K v (r) : r > 0} is a 
modified Bessel function which is positive and smooth in 1Z + , has a pole at the 
origin, and decays exponentially (Abramowitz and Stegun (1970)). Consequently, 
(f c is a non-negative admissible CND1 function. Further, the exponential decay of 
(p c ensures that the symbol function 

*c(0= E ^c(||e + 27rfc||) (4.4.2) 

kez d 
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is continuous for £ e 7£ d \ 2-7rZ d . Therefore, given any finite subset N of Z d , 
Theorem 4.3.10 and Proposition 4.2 imply that the distance matrix A N has every 
eigenvalue bounded away from zero by at least 

»c = \<Pc(\\™ + 2nk\\)\, (4.4.3) 

kez d 

where e = [1, 1, . . ., 1] T G 1Z d . Moreover, Corollary 4.3.11 shows that this bound 
is optimal. 

It follows from (4.4.3) that \x c — > as c — > oo, because of the exponential 
decay of the modified Bessel functions for large argument. For example, in the 
univariate case we have the formula 



He = (4c/tt) \k 1 (ctt) + Ki(3ctt)/3 + Ki(5ctt)/5 



+ • 



and Table 4.1 displays some values of fi c . Of course, a practical implication of 
this result is that we cannot expect accurate direct solution of the interpolation 
equations for even quite modest values of c, at least without using some special 
technique. 



c 


Optimal bound 


1.0 


4.319455 x 10" 2 


2.0 


2.513366 x 10" 3 


3.0 


1.306969 x 10" 4 


4.0 


6.462443 x 10" 6 


5.0 


3.104941 x 10" 7 


10.0 


6.542373 x 10" 14 


15.0 


2.089078 x 10" 20 



Table 4.1: The optimal bound on the smallest eigenvalue as c — > oo 

The optimal bound is achieved only when the numbers of centres is infi- 
nite. Therefore it is interesting to investigate how rapidly |A^ in | converges to the 
optimal lower bound as \N\ increases. Table 4.2 displays |A^ in | = £t c (n), say, for 
the distance matrix (y? c (||j — k\\))™^ for several values of n when c = 1. The 
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third column lists close estimates of /U c (n) obtained using a theorem of Szego (see 
Section 5.2 of Grenander and Szego (1984)). Specifically, Szego's theorem provides 
the approximation 

fi c (n) « <J c (n + n/n), 

where cr c is the function defined in (4.4.2). This theorem of Szego requires the 
fact that the minimum value of the symbol function is attained at n, which is 
inequality (4.3.7). Further, it provides the estimates 

A fc+ i w a c (n + kn/n), fc = 1, . . . , ra - 1, 

for all the negative eigenvalues of the distance matrix. Figure 4.1 displays the 
numbers { — 1/Xk ■ k = 2, ...,n} and their estimates { — l/a(iv + kiv/n) : k = 
l,...,n — 1} in the case when n = 100. We see that the agreement is excel- 
lent. Furthermore, this modification of the classical theory of Toeplitz forms also 
provides an interesting and useful perspective on the construction of efficient pre- 
conditioned for the conjugate gradient solution of the interpolation equations. We 
include no further information on these topics, this last paragraph being presented 
as an aperitif to the paper of Baxter (1992c). 



n 


Hi(n) 




ai(iv + iv/n) 




100 


4.324685 x 10" 


-2 


4.324653 x 10" 


-2 


150 


4.321774 x 10" 


-2 


4.321765 x 10- 


•2 


200 


4.320758 x 10" 


-2 


4.320754 x 10" 


■2 


250 


4.320288 x 10" 


-2 


4.320286 x 10" 


•2 


300 


4.320033 x 10" 


-2 


4.320032 x 10- 


-2 


350 


4.319880 x 10" 


-2 


4.319879 x 10" 


•2 



Table 4.2: Some calculated and estimated values of A^- when c = 1 
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Key: + - calculated, x - estimated 




Figure 4.1. Spectral estimates for a distance matrix of order 100 

4.5. A stability estimate 

The purpose of this last note is to derive an optimal inequality of the form 

2 



TZ d 



2^ vm\\n-3\ 

jez d 



jez d 



where {yj)j e z d is a rea l sequence of finite support such that Yljez d Vi = 0> an< ^ 
ip: [0, oo) 1Z belongs to a certain class of functions including the multiquadric. 
Specifically, this is the class of admissible CND1 functions. These functions have 
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generalized Fourier transforms given by 

/•oo 

m\\) = - exp(-||e|| 2 /*)d/*(t), 
Jo 



(4.5.1), 



where dfi is a positive (but not finite) Borel measure on [0, oo). A derivation of 
this expression may be found in Theorem 4.2.6 above. 

Lemma 4.5.1. Let {yj)j^z d be a zero- summing sequence and let <p: [0, oo) — > 1Z 
be an admissible CND1 function. Then we have the equation 



Yl yM\\ x -j\ 

jez d 



dx = (2n)~ d [ | J2 VjeMijOlMOdZ, 
J[0W jezd 



where a(0 = E keZ d \m + ^k\\)\ 2 - 

Proof. Applying the Parseval theorem and dissecting 1Z d into copies of the cube 
[0, 27r] d , we obtain the equations 

/ | yMW* - j\\)\ 2 dx 

Jnd je z d 

= (2n)- d I \Y.V3^Mm\ 2 \m\\)\ 2 di 

= W d I I E I/iex P (yOri^(||e + 27rfc||)| 2 de 

kezd J[0,2*]<* jeZd 

= (2n)~ d [ | J2 VjeMijOlMOdd, 

J[0,2ir] d d 



jez d 

where the interchange of summation and integration is justified by Fubini's theo- 
rem. ■ 

If c(0 — m f° r almost every point £ in [0, 2ir] d , then the import of Lemma 
4.5.1 is the bound 



K d 



E yM\\ x ~j\ 

jez d 



dx > m Uj. 

jez d 



63 



Norm estimates for Toeplitz distance matrices I 

We shall prove that we can take m = cr(7re), where e = [1, 1, . . . , 1] T e 1Z d . Further, 
we shall show that the inequality is optimal if the function a is continuous at the 
point Tie. 

Equation (4.5.1) is the key to this analysis, just as before. We see that 



pOO pOO 

\m\\)\ 2 = / / exp(-^ii 2 (tr 1 +t 2 - 1 ))«i)«2), 

Jo Jo 

whence, 

poo poo 

a & = / / E ex PHI£ + 2nkf(t^ + t^ 1 )) dfifa) dn(t 2 ), (4.5.2), 
Jo Jo kez d 

where the interchange of summation and integration is justified by Fubini's theo- 
rem. 

Now it is proved in Lemma 4.1.8 that 

J2 exp(-A||£ + 27rfc|| 2 ) > exp(-A||7re + 27rfc|| 2 ), 
kez d kez d 

for any positive constant A. Therefore equation (4.5.2) provides the inequality 

poo poo 

<r(£)> / / E exp(-||7re + 27r/c|| 2 (^ 1 +t2 1 ))^(^i)^(^), 
Jo Jo kez d 

= o-(ne), 

which is the promised value of the lower bound m on a mentioned above. Thus 
we have proved the following theorem. 

Theorem 4.5.2. Let (yj)j EZ d, <p and a be as defined in Lemma 1. Then we have 
the inequality 



vM\\ x -j\ 

jez d 



dx > a (Tie) ^2 y]. 

jez d 



The proof that this bound is optimal uses the technique of Theorem 4.2.11. 

Theorem 4.5.3. The inequality of Theorem 2 is optimal if a is continuous at ne. 

64 



Norm estimates for Toeplitz distance matrices I 

Proof. The condition that ip be admissible requires the existence of the limit 
lim||^|| ll£l| d+ V(ll£ll)-Let m be a positive integer such that 2m > d + 1 and 
let us define a sequence {(Uj )j£Z d '■ n = 1,2, . . .} by 



jez d 



2m 



cT 1 ]Tsin 2 (^/2) K n (e-7re), 



us to 



where K n denotes the multivariate Fejer kernel. The standard properties of the 
Fejer kernel needed for this proof are described in Lemma 4.1.2. They allow 
deduce that (y^)j^z d is a zero-summing for every n. Further, we see that 

£ \yf\ 2 = (2n)~ d [ K^-neXd-^^&m) 2 ™^ 

£2* pi 

= 1, for n > 4m. 
Finally, m has been chosen so that the function 

(\ 2m 
d-'p^^m) a(0:£e[O,27r] d } 

is continuous. Therefore, we have 



lim 



yM\\ x ~j\ 

jez d 



dx 



= lim (2n)~ d / ^(^-Tre)^" 1 Vsin 2 (^/2)) 2m a(0^ 

n ^°° J[0,2n] d pi 

= a(ne), 

using the fact that a is continuous at Tie and standard properties of the Fejer 
kernel. ■ 

4.6. Scaling the infinite grid 

Here we consider the behaviour of the norm estimate given above when we scale 
the infinite regular grid. 
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Proposition 4.6.1. Let r be a positive number and let (dj)j eZ d be a real sequence 
of finite support. Then 



a jak eM-\\rj-rk\\ 2 ) = (27r)" d / I ]T a.expz^ V^df, (4-6.1 







j,kez d - L ">-»j i6 ^ 



= e- ||rfc " 2 expz/cC £ G (4.6.2) 



Proof. Section 4.2 provides the equation 

^2 aja k exp(-\\rj - rkf) 

j,kez d 



= (2n)~ d [ I "jexpitt \^/r 2 ) d/2 £ exp(-|K + 27rfc|| 2 /4r 2 ) 

J[0,2^ jeZd keZd 

(4.6.3) 

Further, the Poisson summation formula gives the relation 

(2n) d J2 e M~U + 27rfc|| 2 /4r 2 ) = (4nr 2 ) d / 2 expz/c£. (4.6.4) 

kez d kez d 

Substituting (4.6.4) into (4.6.3) yields equations (4.6.1) and (4.6.2). | 

The functions E^ and Er are related in a simple way. 
Lemma 4.6.2. We have the expression 

4 d) (0=U E r 1) ^). (4.6.5) 
fc=i 

Proof. This is a straightforward consequence of (4.6.2). | 
Applying the theta function formulae of Section 4.2 yields the following 

result. 

Lemma 4.6.3. 

oo 

= J](l - e- 2fcr2 )(l + 2e"( 2fc - 1 ^ 2 cos£ + e^ 4 *" 2 ^). (4.6.6) 
fc=i 
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Proof. The Theta function #3 of Jacobi type is given by 



00 

L.2 



83(2, q) = 1 + 2 cos2kz, q,z G C, \q\ < 1, 

oo (4-6-7) 
= - ? 2fc )(l + 2Q 2fc " 1 cos 2,2 + Q 4fc " 2 ), 

fc=i 

2 

which equations are discussed in greater detail in Section 4.2. Setting q = e~ r we 
have the expressions 

00 

EP(Z) = 6 3 (Z/2,q) = JJ(l-e- 2fcr2 )(l + 2e- (2fc - 1)r2 cos^ + e- (4fc - 2)r2 ). (4.6.8) 

k=i 

The proof is complete. | 
Now Section 4.3 provides the inequality 

£ «^fcexp(-||rj - r£;|| 2 ) > E^^e) £ a 2 , (4.6.9) 

j,kez d jez d 

where e = [1, . . . , 1] T G 7?. d . Using equation (4.6.6), we see that 

= JJ(1 - e- 2fcra )(l - 2e-( 2fc " 1 ^ 2 + e-( 4fc " 2 ^ 2 ) 
fc=i 

= J](l-e- 2fc7 ' 2 )(l-e-( 2fc - 1 ^ 2 ) 2 , 



(4.6.10) 



fc=i 



which implies that {E^ : r > 0} is an increasing function. Further, it is a conse- 
quence of (4.6.5) that {Er (7re) : r > 0} is also an increasing function. We state 
these results formally. 

Theorem 4.6.4. Let r > s > 0. Then we have the inequality 

inf aja fc exp( — \\rj — rk\\ 2 ) > inf ajafc exp( — — s/c|| 2 ), (4.6.11) 

j,kez d j,kez d 

where the infima are taken over the set of real sequences of finite support. 
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In fact we extend the given analysis to a class of functions including the 
multiquadric. The appropriate definitions and theorems form Section 4.3, but the 
key result is Theorem 4.3.6: Under suitable conditions, the function ip: [0, oo) — > TZ 
possesses the generalized Fourier transform 

POO 

m\\) = - exp(-||^|| 2 /4t)t- 1 dp(t), (4.6.11) 
Jo 

where 

/•oo 

<p(r) = <p(0)+ (l-e^ 2 *)* -1 <*//(*)> (4.6.12) 
Jo 

and \x is a positive Borel measure such that f Q dji{t) < oo and t _1 d[i(t) < oo. 
Now the function ip r : x t— > ip(\\rx\\) has the Fourier transform 

MU\\) = m\\Mr- d . (4.6.13) 

Further, the associated symbol function is defined by the equation 

M0= E \MU + 2nk\\)\, (4.6.14) 

kez d 

and so (4.6.13) implies the expression 

/•OO 

ff r(0 = / r~ d V exp(-||£ + 2nkf/4tr 2 )(n/t) d / 2 t- 1 d/*(f). (4.6.15) 

Jo kez d 

Using the Poisson summation formula, we have 

(2n) d exp(-||£ + 27r£;||74tr 2 ) = (4tr 2 n) d / 2 e'^^e^. (4.6.16) 
kez d kez d 

Consequently we have 

\n/t) d ' 2 eM-U + 2nk\\ 2 /4tr 2 )= £ e ~^\^ = E% 2 (0, 
kez d kez d 

(4.6.18) 



r -d, 



providing the equation 

/•oo 

a r (ne)= / E%> 1/3 (ire)^ 1 d»(t), (4.6.18) 
Jo 

and so {a r (7re) : r > 0} is an increasing function. 
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Appendix 

I do not like stating integral representations such as Theorem 4.3.5 without in- 
cluding some explicit examples. Therefore this appendix calculates da for ip{r) = r 
and (p(r) = (r 2 + C 2 ) 1 / 2 , where c is positive and we are using the notation of 4.3.5. 
For (p(r) = r the key integral is 

1 f°° 

r(~) = J {e~ u - 1) du, (Al) 

which is derived in Whittaker and Watson (1927), Section 12.21. Making the sub- 
stitution u = r 2 t in (Al) and using the equations 7T 1 / 2 = T(l/2) = — r(— 1/2)/2 
we have 

(e-^-l)t- 3 / 2 ^, 

that is 

r = J (l-e'^^f- 1 (4nt)- 1/2 dt. (A2) 

Thus the Borel measure is da\(t) = (47rf) -1 / 2 dt and f Q da\{t) = t~ x da\(t) = 
n 1 / 2 . 

The representation for the multiquadric is an easy consequence of (A2). 
Substituting (r 2 + c 2 ) 1 / 2 and c for r in (A2) we obtain 

(r 2 + c2) l/2 = J ^ _ e -(r 2 +c 2 )t j f -l ( 47rt )-l/2 dt (A3) 

and 

POO 

c = J (l-e-^jt- 1 (Ant)' 1 ' 2 dt, (A4) 
respectively. Subtracting (A4) from (A3) provides the formula 

(l - e~ rH ) t- 1 e- c2t (4yrt)- 1 /2 dt . (A5) 
Hence the measure is da2(t) = e~ c2t (4nt)~ 1 / 2 dt. 
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5.1. Introduction 

Let <p: lZ d — > 1Z be an even continuous function of at most polynomial growth. As- 
sociated with this function is a symmetric bi-infinite multivariate Toeplitz matrix 

*=fa(j-*)W-- (5.1.1) 
Every finite subset I = (ij)™=i of Z d determines a finite submatrix of $ given by 

d> 7 := (pfe - fc fc))? fc=1 . (5.1.2) 

We are interested in upper bounds on the £ 2 -norm of the inverse matrix that 
is the quantity 

||$7 1 || = l/min{||x|| 2 : \\&ix\\ 2 = 1, (5.1.3) 

where ||x||2 = X^e/^j for £ = ( x j)jei- The type of bound we seek follows the 
pattern of results in the previous chapter. Specifically, we let (p be the distributional 
Fourier transform of ip in the sense of Schwartz (1966), which we assume to be a 
measurable function on 1Z d . We let e := (1, . . . , 1) T G 1Z d and set 

r := \^e + 2nj)\ (5.1.4) 

jez d 

whenever the right hand side of this equation is meaningful. Then, for a certain 
class of radially symmetric functions, we proved in Chapter 4 that 

W^W < 1/t# (5.1.5) 

for every finite subset / of Z d . Here we extend this bound to a wider class of 
functions which need not be radially symmetric. For instance, we show that (5.1.5) 
holds for the class of functions 

ip(x) = (||x||i +c) 7 , xeTZ d , 
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where ||x||i = X^=i \ x j\ ^ s the ^i-norm of x, c is non-negative, and < 7 < 1. 

Our analysis develops the methods of Chapter 4. However, here we empha- 
size the importance of certain properties of Polya frequency functions and Polya 
frequency sequences (due to I. J. Schoenberg) in order to obtain estimates like 
(5.1.5). 

In Section 2 we consider Fourier transform techniques which we need to 
prove our bound. Further, the results of this section improve on the treatment 
of the last chapter, in that the condition of admissibility (see Definition 5.3.2) is 
shown to be unnecessary. Section 3 contains a discussion of the class of functions 
if for which we will prove the bound (5.1.4). The final section contains the proof 
of our main result. 



5.2. Preliminary facts 

We begin with a rather general framework. Let (p: 1Z d — > 1Z be a continuous func- 
tion of polynomial growth. Thus p> possesses a distributional Fourier transform 
in the sense of Schwartz (1966). We shall assume (p is almost everywhere equal 
to a Lebesgue measurable function on 7Z d , that is we assume (p to be the sum 
of a measurable function and a tempered distribution whose support is a set of 
Lebesgue measure zero. Given a nonzero real sequence (yj)j eZ d of finite support 
and points {x^)j eZ d in 1Z d , we introduce the function F:1Z d —¥ 1Z given by 

F(x)= Vjyk<p{x + x j ~x k ), xell d . (5.2.1) 

j,kez d 

Thus 

^(0)= Yl VjyMx j ~x k ), (5.2.2) 

j,kez d 

which is the quadratic form whose study is the object of much of this dissertation. 
We observe that the Fourier transform of F is the tempered distribution 

^(0=|E W etei€ |V(0, Z£K d . (5.2.3) 

jez d 
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Further, if F is an absolutely integrable function, then we have the equation 

F(0) = (2n)- d [ F(Odf, (5.2.4) 
Jn d 

since F is the inverse distributional Fourier transform of F and this coincides with 
the classical inverse transform when F e L 1 {TZ d ). In other words, we have the 
equation 

y 3 y k ^ -x k ) = {2n)~ d j^Vje^^mdi (5.2.5) 

j,kez d nd jez d 

when F is absolutely integrable. If we make the further assumption that (p is one- 
signed almost everywhere on lZ d : and the points {x^)j eZ d form a subset of the 
integers Z d , then it is possible to improve (5.2.5). First observe that dissecting 1Z d 
into 2-7r-integer translates of the cube [0,2-7r] d provides the relations 



kez d nd jez d 

= [ I E Vi* 1 '* ^ + 27rfc ) * (5.2.6) 



keZ d J[0,2n]- jeZ d 

(27T)"" 



J[0,27T]dl d I 



jez d 
where 

<K0 = E £(£ + 2 **). £eft d , (5-2.7) 

kez d 

and the monotone convergence theorem justifies the exchange of summation and 
integration. Further, we see that another consequence of the condition that (p be 
one-signed is the bound 

|E%' e ^| 2 |^)l<oo 
jez d 

for almost every point £ G [0,27r] d , because the left hand side of (5.2.6) is a 
fortiori finite. This implies that a is almost everywhere finite, since the set of all 
zeros of a nonzero trigonometric polynomial has measure zero. This last result 
is well-known, but we include its short proof for completeness. Following Rudin 
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(1973), we shall say that a continuous function f:C d ^rC is an entire junction of d 
complex variables if, for every point (w\, . . . , w d ) € C d and for every j E {1, . . . , d}, 
the mapping 

C 3 z ^ f(wi, . . . , Wj-i, z, Wj+i, ...,w d ) 
is an entire function of one complex variable. 

Lemma 5.2.1. Given complex numbers (a_j)?=i an< ^ a se ^ of distinct points (x 3 )™ 
in lZ d , we let p: lZ d — >■ C be the function 

n 
3 = 1 

Then p enjoys the following properties: 

(i) p is identically zero if and only if aj =0, 1 < j < n. 
(ii) p is nonzero almost everywhere unless aj = 0, 1 < j < n. 

Proof. 

(i) Suppose p is identically zero. Choose any j e {1, . . . , n} and let fy. lZ d — > 1Z be 
a continuous function of compact support such that fj(x k ) = djk for 1 < k < n. 
Then 

n „ n 

«i = Va fc /,(a; fc ) = (27r)- d / Va/^(()^ = 0. 
The converse is obvious. 

(ii) Let /: C d — > C be an entire function and let 

Z = {x G 1Z d : f(x) = 0}. 

If voldZ is a set of positive Lebesgue measure in 1Z d , then we shall prove that / is 
identically zero, which implies the required result. 

We proceed by induction on the dimension d. If d = 1 and voliZ > 0, 
then / is an entire function of one complex variable with uncountably many zeros. 
Such a function must vanish everywhere, because every uncountable subset of C 
possesses a limit point. Now suppose that the result is true for d — 1 for some 
d > 2. Fubini's theorem provides the relation 

< vol d Z = / vol 1 Z(x 2 , ■ ■ ■ , x d ) dx 2 ■ ■ ■ dx d , 
Jn d - 1 
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where 

Z(x 2 , ...,x d ) = {x 1 eTZ:(x 1 ,...,x d )e Z}. 

Thus there is a set, X say, in 7^. d_1 of positive (d — 1) -dimensional Lebesgue 
measure such that vol\Z(x2, • • • , x d ) is positive for every (x 2 , • • • , Xd) G X, and 
therefore the entire function C 3 z i— > f(z,X2,---,Xd) vanishes for all z G C, 
because Z(x 2 , ■ ■ ■ , Xd) is an uncountable set. Thus, choosing any z\ G C, we see 
that the entire function of d — 1 complex variables defined by 

(z 2 , ...,Zd) >-> f(zi, Z2,..., Z d ), (Z2, ...,Zd) G C d_1 , 

vanishes for all (z2, ■ ■ ■ , z d ) in X, which is a set of positive (d — 1)- dimensional 
Lebesgue measure. By induction hypothesis, we deduce that 

f(z 1 , z 2 , . . ■ , z d ) = for all z 2 , . . . , z d G C, 

and since z\ can be any complex number, we conclude that / is identically zero. 
Therefore the lemma is true. ■ 

We can now derive our first bounds on the quadratic form (5.2.2). For 
any measurable function g: [0,27r] d — > 1Z, we recall the definitions of the essential 
supremum 

ess sup g = inf{c G 1Z : g(x) < c for almost every x G [0, 27r] d } (5.2.8) 
and the essential infimum 

ess inf g = sup{c G 1Z : g(x) > c for almost every x G [0, 2n] d }. (5.2.9) 
Thus (5.2.6) and the Parseval relation provide the inequalities 

ess inf a ^ y| < ^ y^y^ix* — x k ) < ess sup a ^ y|. (5.2.10) 

jez d j,kez d jez d 

Let V be the vector space of real sequences (yj)j EZ d of finite support for which 
the function F of (5.2.3) is absolutely integrable. We have seen that (5.2.10) is 
valid for every element (yj)j e z d °f V. Of course, at this stage there is no guarantee 
that V 7^ {0} or that the bounds are finite. Nevertheless, we identify below a case 
when the bounds (5.2.10) cannot be improved. This will be of relevance later. 

74 



Norm estimates for Toeplitz distance matrices II 

Proposition 5.2.2. Let P be a nonzero trigonometric polynomial such that the 
principal ideal I generated by P , that is the set 



X = {PT : T a real trigonometric polynomial} , (5.2.11) 

consists of trigonometric polynomials whose Fourier coefficient sequences are ele- 
ments ofV. Further, suppose that there is a point rj at which a is continuous and 
P(rj) 7^ 0. Then we can find a sequence {{y^)jez d '■ n = 1,2, . . .} in V such that 

lim E yfvP^j-k)/j2^ 2 = ^)- ( 5 - 2 - 12 ) 

j,kez d jez d 

Proof. We recal Section 4 and recall that the nth degree tensor product Fejer 
kernel is defined by 

■= II ^-fi'l = ^ E ^ 2 = : WOf, t e ^ (5-2-13) 



J — l 3 -" k<EZ a 

0<k<en 



where e = (1, . . . , 1) T G ft d and L n (f ) = n~ d / 2 Eo<fc<en Then the function 
P(-)L n (-—rj) is a member of X and we choose (y^)j eZ d to be its Fourier coefficient 
sequence. The Parseval relation provides the equation 

E [yf? = W d [ P 2 (0K n (t - V) dt (5.2.14) 

jeZ d J[0,2n] d 

and the approximate identity property of the Fejer kernel (Zygmund (1988), p. 86) 
implies that 

P 2 ( V ) = lim (27r)" d / P 2 (0K n (H - rj) 

_ . \ J 5.2.15 
= lim E [1/ n) ] 2 - 

Further, because a is continuous at r\, we also have the relations 



™ ^ ^ (5.2.16) 
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the last line being a consequence of (5.2.6). Hence (5.2.15) and (5.2.16) provide 
equation (5.2.12). ■ 

Corollary 5.2.3. If a attains its essential infimum (resp. supremum) at a point 
of continuity, and if we can find a trigonometric polynomial P satisfying the con- 
ditions of Proposition 5.2.2, then the lower (resp. upper) bound of (5.2.10) cannot 
be improved. 

Proof. This is an obvious consequence of Proposition 5.2.2. ■ 

We now specialize this general setting to the following case. 

Definition 5.2.4. Let G: lZ d — > 1Z be a continuous absolutely integrable function 
such that G(0) = 1 for which the Fourier transform is non-negative and absolutely 
integrable. Further, we require that there exist non-negative constants C and k 
for which 

|1 -G(x)\ < C\\x\\ K , xeTZ d . (5.2.17) 

We let Q denote the class of all such functions G. 

Clearly the Gaussian G(x) = exp( — ||x|| 2 ) provides an example of such a 
function. The next lemma mentions some salient properties of Q which do not, 
however, require (5.2.17). 

Lemma 5.2.5. Let G G Q . 

(i) G is a symmetric function, that is 

G(x) = G(-x), xeIZ d . (5.2.18) 

(ii) 

\G(x)\ < 1, xeTZ d . (5.2.19) 

(Hi) G is a positive definite function in the sense of Bochner. In other words, for 
any real sequence (yj)j EZ d of finite support, and for any points {x^)j eZ d in 
1Z d , we have 

VjVkGix^ - x k ) > 0. (5.2.20) 

j,kez d 
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Proof. 

(i) Since G is real-valued we have 



2i G(x)smxtdx = G(Z)-G(-Z)£K, £ G K d , 
Jn d 

which is a contradiction unless both sides vanish. Thus G is a symmetric 
function. However, G must inherit this symmetry, by the Fourier inversion 
theorem. 

(ii) The non-negativity of G provides the relations 

\G(x)\= (2n)~ d [ G(Oe- te *de <(27r)- d / G(£) df = G(0) = 1. 

(hi) The condition G G L 1 (TZ d ) implies the validity of (5.2.5) for </> replaced by 
G, whence 

, , y. 2 



yj y k G(xi - x k ) = (2n)~ d [\Y1 Vi 

j,kez d nd jez d 



G(0 dH > 0, 



as required. 



We remark that the first two parts of Lemma 5.2.5 are usually deduced from 
the requirement that G be a positive definite function in the Bochner sense (see 
Katznelson (1976), p. 137). We have presented our material in this order because 
it is the non- negativity condition on G which forms our starting point. 

Given any G G Q, we define the set A(G) of functions of the form 

poo 

(p(x) = c+ [1 - G(t 1/2 x)]t _1 da(t), xe1l d , (5.2.21) 
Jo 

where c is a constant and a: [0, oo) — > 1Z is a non-decreasing function such that 

/OO />1 
t~ x da(t) < oo and J t K/2 ~ x da(t) < oo. (5.2.22) 

Let us show that (5.2.21) is well-defined. Inequality (5.2.19) implies the 



bound 



/oo />oo 
1 - G(t 1/2 x) t' 1 da(t) < 2 J r 1 da(t) < oo. 
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Moreover, applying condition (5.2.17) we obtain 

/ l-G(t 1/2 x) t^dctit) < C\\x\\ R [ W 2 - 1 dct(t) < oo. (5.2.23) 
Jo Jo 

Therefore the integral of (5.2.21) is finite and ip is a function of polynomial growth. 

A simple application of the dominated convergence theorem reveals that ip is also 

continuous, so that we may view it as a tempered distribution. 

The following definition is convenient. 

Definition 5.2.6. We shall say that a real sequence (yj)j e z d of finite support is 
zero- summing if J2jez d Vj = 0- 

An important property of A(G) is that it consists of conditionally negative 
definite functions of order 1 on 1Z d , that is whenever <p G A(G) 

V3VW{x j ~x k )<0 (5.2.24) 

j,kez d 

for every zero-summing sequence (yj)j eZ d an d for any points (xj)j eZ d m l - 
Indeed, (5.2.21) provides the equation 

/•OO 

Y, ViVk^-x k ) = - / J2 VjVkGit 1 ' 2 ^ -x^t-Uait), (5.2.25) 
j,kez d Jo j,kez d 

and the right hand side is non-positive because G is positive definite in the Bochner 
sense (Lemma 5.2.5 (iii)). 

We now fix attention on a particular element G G Q and a function <p G 

A(G). 

Theorem 5.2.7. Let (yj)j e z d be a zero-summing sequence that is not identically 
zero. Then, for any points (x J )j €2 d in lZ d , we have the equation 

yj yMxi-x k ) = -(2n)~ d f J V^^HiOdt;. (5.2.26) 

j,kez d Jnd jez d 

where 

POO 

H(0= G(^/t 1/2 )t- d/2 - 1 rfa(t), ieTZ d . (5.2.27) 
Jo 

Furthermore, this latter integral is finite for almost every £ G lZ d . 
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Proof. Applying the Fourier inversion theorem to G in (5.2.25), we obtain 



j,kez d 

"OO 



/>oo r 2 

= -(2n)- d / V yj exp(it 1/2 r]x j ) Girfir 1 dr] dot(t) 

Jo Jni\ jezd 

= -{2-K)~ d / \Y,Vi^ G{i/t^)r d ' 2 - l dida{t), 

Jo Jn d 1 ,^~ d 



(5.2.28) 



jez d 

where we have used the substitution £ = t 1 / 2 ??. Because the integrand in the last 
line is non-negative, we can exchange the order of integration to obtain (5.2.26). 
Of course the left hand side of (5.2.26) is finite, which implies that the integrand of 
(5.2.26) is an absolutely integrable function, and hence finite almost everywhere. 
But, by Lemma 5.2.1, | J2j Uje^^l 2 ^ for almost every £ e 1Z d if the sequence 
(Uj)jez d is non-zero. Therefore H is finite almost everywhere. ■ 

Corollary 5.2.8. The hypotheses of Theorem 5.2.7 imply the equation 

F(x) = -(2n)~ d [ I V y^* V (£)e^ df, x G TZ d , (5.2.29) 
Jnd jez* 

where F is given by (5.2.1). Consequently, = —H{£) for almost every £ G lZ d , 
that is 

/•OO 

m = - / G(f/t 1/2 )t -d/2-1 da(f). (5.2.30) 
Jo 

Proof. It is straightforward to deduce the relation 

= -(2yr)- d / / £ e^Git/t 1 ' 2 )^' 2 - 1 d£da(t), 

Jo Jn dl jezd 

which is analogous to (5.2.28). Now the absolute value of this integrand is precisely 
the integrand in the second line of (5.2.28). Thus we may apply Fubini's theorem 
to exchange the order of integration, obtaining (5.2.29). 

Next, we prove that £ i-> — | ^\ yje lx3 ^\ 2 H(^) is the Fourier transform of 
F. Indeed, let ift: lZ d — > 1Z be any smooth function whose partial derivatives enjoy 
supra-algebraic decay. It is sufficient (see Rudin (1973)) to show that 

f $(x)F(x) dx=- f V(0| E Vi e ^ ^ d ^ 

Jn d Jn d 1 ~z d 

jez d 
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Applying (5.2.29) and Fubini's theorem, we get 

f ^{x)F{x)dx = -{2>K)- d [ [ ^(x)\Y,yje lx ^ 2 e^H(OdCdx 



f I E Vi^ 2 H(0 (W" I tWe** dx) d£ 

I E ^"'" ^i^ciodt. 



jez d 

which establishes (5.2.30). However, we already know that the Fourier transform 
F(£) is almost everywhere equal to | J2j yj elxJ ^\ 2 0(O- By Lemma 5.2.1, we know 
that J2jVj elxJ ^ 7^ f° r almost all £ G 1Z d , which implies that tp = —H almost 
everywhere. ■ 

5.3. Polya frequency functions 

For every real sequence {cij)JL 1 and any non-negative constant 7 such that < 



00 

2 



E(z) = e-^ z Y[(l-a 2 jZ 2 ), zeC. (5.3.1) 

This is an entire function which is nonzero in the vertical strip 

\$lz\ < p := 1/ sup{|aj I : j = 1,2,...}. 

It can be shown (Karlin (1968), Chapter 5) that there exists a continuous function 
Allien such that 

/ A(t)e~ zt dt = — \Uz\<p. (5.3.2) 
Jn E \ z ) 

This function A is what Schoenberg (1951) calls a Polya frequency function. We 
have restricted ourselves to functions A which are even, that is 



A(£) = A(-£), tell. 
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Also, E(0) = 1 implies that 



L 



A(t)dt = l. (5.3.4) 



According to (5.3.1) the Fourier transform of A is given by 



c2 



We see that A(-)/A(0) is a member of the set Q described in Definition 5.2.4 for 
d = 1, and therefore Lemma 5.2.5 is applicable. In particular, 

|A(t)| < A(0), tell. (5.3.6) 

However, much more than (5.3.6) is true. Schoenberg (1951) proved that 

detCA^-yfc^i^O (5-3.7) 

whenever x\ < • ■ • < x n and yi < • • • < y n - This fact will be used in an essential 
way in Section 5.4. For the moment we use it to improve (5.3.6) to 

A(£) e [0, A(0)], teTZ. (5.3.8) 

Let V denote the class of functions k-.n — > n that satisfy (5.3.2) for some 7 > 
and sequence (aj)J^ 1 satisfying < 7 + J2jLi a j < 00 • For any positive a the 
function 

Sa(t) = 7 ±-e-W a \, ten, (5.3.9) 

Zi\CL\ 

is in V since 

/ S a (t)e- Zt dt = l -^r, \Uz\ < 1/a. (5.3.10) 

Jn 1 - a z 

Let 8 = {S a : a > 0}. These are the only elements of V that are not in C 2 {n), 
because all other members of V have the property that kit) = 0(t~ 4 ) as \t\ — > 00. 
Hence there exists a constant k such that 

|A(0) - A(£)| < nt 2 , for t e n and A e V\E, (5.3.11) 

or 

|A(0) — A(*)| < k|*|, ten, keS. (5.3.12) 
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We note also that every element of V decays exponentially for large argument (see 
Karlin (1968), p. 332). 

We are now ready to define the multivariate class of functions which interest 
us. Choose any A 1 ,...,A (J eP and define 

G W = lll^' x = ( Xl ,...,x d )en d . (5.3.13) 

According to (5.3.11) and (5.3.12), there is a constant C > such that 

1 - G(x) < C\\x\\l, (5.3.14) 

when Aj ^ £ for every factor Aj in (5.3.11). However, if Aj E £ for every j, then 
we only have 

1 -G(x) < C\\x\\ 2 , (5.3.15) 

for some constant C. We are unable to study the general behaviour at this time. 
Remarking that the Fourier transform of G is given by 

d ^ = II x|f £= •••,&) e n d , (5.3.16) 

we conclude that G is a member of the class Q of Definition 5.2.4. Moreover, we can 
now construct the set A(G). To this end, let a: [0, oo) — > 71 be a non-decreasing 
function such that 

/oo 
t _1 rfa(t)<oo, (5.3.17) 

and for any constant c E TZ define (p: [0, oo) — > 7£ by (5.2.21). Thus we see that as 
long as we require the measure da to satisfy the extra condition 

/ t" 1/2 da(t) < oo (5.3.18) 

whenever one of the factors in (5.3.11) is an element of £, then ip is a continuous 
function of polynomial growth and the results of Section 2 apply. We let C denote 
the class of all such functions, for all G E Q. 
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Let us note that C contains the following important subclass of functions. 
In 1938, I. J. Schoenberg proved that a continuous radially symmetric function 
ip: lZ d — > 1Z is conditionally negative definite of order 1 on every lZ d if and only if 
it has the form 

(p(x) = (p(0) + J (l -exp(-t||x|| 2 ))t _1 rfa(t), x G 1Z d , 

where a: [0, oo) — > 1Z is a non-decreasing function satisfying (5.3.17). In this case, 
the Gaussian is clearly of the form (5.3.13), implying that we do indeed have a 
subclass of C. Thus we have established Theorem 5.2.7 and Corollary 5.2.8 under 
weaker conditions than those assumed in Chapter 4. 
Our class C also contains functions of the form 

V(x) = c + J^ (l -exp(-t 1/2 ||x||i))t" 1 rfa(t), x G 7Z d , 

where a: [0, oo) — > 1Z is a non-decreasing function satisfying (5.3.17) and (5.3.18), 
and ||x||i = ^7=i \ x j\ f° r x = { x ii • • • i x d) e F° r instance, using the easily 
verified formula 

which is valid for 5 > and —1/2 < 7 < 0, we see that (p(x) = (5 + ||x||i) r , for 
5 > and < r < 1, is in our class C. 

Although it is not central to our interests in this section, we will discuss 
some additional properties of the Fourier transform of a function (p G C. First, 
observe that (5.3.5) implies that A is a decreasing function on [0, 00) for every A 
in V. Consequently every G G Q satisfies the inequality < G(rf) for £ > 77 > 0. 
This property is inherited by the function H of (5.2.27), that is 

H(£) < H{rj) whenever £ > 77 > 0, (5.3.19) 

which allows us to strengthen Theorem 5.2.7. 
Proposition 5.3.1. H is continuous on (1Z \ {0}) d . 
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Proof. We first show that H is finite on (1Z\ {0}) d . We already know that (p = —H 
almost everywhere, which implies that every set of positive measure contains a 
point at which H is finite. In particular, let 8 be a positive number and set Us = 
{£ G K d : < f j < 5, j = 1, . . . , d}. Thus there is a point r] e Us such that 
H(rj) < oo. Applying (5.3.19) and recalling that if is a symmetric function, we 
deduce the inequality 

H(£) < H( V ) < oo, £ G F 5 , (5.3.20) 

where i 7 ^ := {£ G 7£ d : > S, j = 1, . . . , d}. Since 6 > is arbitrary, we see 
that H is finite in (K \ {0}) d . 

To prove that H is continuous in Fs, let (£ n )^Li be a convergent sequence 
in Fs with limit £oo. By (5.3.20), the functions 

{t^G(t n t- 1/2 )t- d/2 - 1 :n=l,2,...} 

are absolutely integrable on [0, oo) with respect to the measure da. Moreover, they 
are dominated by the cfct-integrable function t >->■ G(i]t~ 1/ ' 2 )t~ d / 2 ~ 1 . Finally, the 
continuity of G provides the equation 

hm G^nt- 1 / 2 )^ 2 - 1 = G(£oo*- 1/2 r d/2 -\ t G [0, oo), 

n— >oo 

and thus lim n ^ OG H(^ n ) = H(^oo) by the dominated convergence theorem. Since 
6 was an arbitrary positive number, we conclude that H is continuous on (71 \ 

{o}) d . m 

The remainder of this section requires a distinction of cases. The first case 
(Case I) is the nicest. This occurs when every factor Aj in (5.3.13) has a posi- 
tive exponent 7j in the Fourier transform formula (5.3.5). We let Case II denote 
the contrary case. Our investigation of Case II is not yet complete, so we shall 
concentrate on Case I for the remainder of this section. 

For Case I we have the bound 

G(0 < e-( 7l ^ + '" +7d ^), ieTZ d , 
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which implies the limit 

limGre*- 1 / 2 )*-^ 2 - 1 = 0, £t^0. 

Thus the function t >->■ G(Z t t~ 1 / 2 )t~ d / 2 ~ 1 is continuous for t e [0, 00) when £ is 
nonzero, which implies that 

f 1 6(&- 1 / 2 )t- d / 2 - 1 da{t) < 00, £^0. 
io 

Moreover, since 

/oo />oo 
G(£,t~ 1/2 )t~ d/2 ~ 1 da(t) < J t" 1 cia(£) < 00, 

we have < 00 for every £ e 7?- d \ {0}. Finally, a simple extension of the proof 
of Proposition 5.3.1 shows that H is continuous on lZ d \ {0}. 

In fact, we can prove that for H e C°°(1Z d \ {0}) in Case I. We observe 
that it is sufficient to show that every derivative of G(^t~ 1 / 2 )t~ d / 2 ~ 1 with respect 
to £ is an absolutely integrable function with respect to the measure da on [0, 00), 
because then we are justified in differentiating under the integral sign. Next, the 
form of G implies that we only need to show that every derivative of A, where A is 
given by (5.3.13) and 7 > 0, enjoys faster than algebraic decay for large argument. 
To this end we claim that for every C < p := 1/ sup{|a J | : j = 1,2,...} there is a 
constant D such that 

Mt + iv) <De~ ie , tell, \r}\<C. (5.3.21) 

To verify the claim, observe that when \q\ < C < |£| we have the inequalities 



< e^^e'^ 2 and |1 + a 2 (£ + ii]) 2 \ > 1 + a 2 (£ 2 - if) > 1. 



Thus, setting M = max{|A(£ + ir])\e^ : |f| < C,\r]\ < C}, we conclude that 
D := max{M, e c 7 } is suitable in (5.3.21). Finally, we apply the Cauchy integral 
formula to estimate the kth derivative. We have 
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where V : [0, 2n] — > C is given by T(t) = re lt and r < C is a constant. Consequently 
we have the bound 



A (fc) (0 < (D/a*)e- 7Ddn {< € - r > 2 '« +r > 3 >, ^ ElZ 



and the desired supra-algebraic decay is established. We now state this formally. 

Proposition 5.3.2. In Case I, the function H of (5.2.27) is smooth for nonzero 
argument. 

Next, to identify —H with (p on 1Z d \ {0} in Case I, we let ip: 1Z d — > 1Z be a 
smooth function whose support is a compact subset of 1Z d \ {0}. By definition we 
have 



(<p,ip)= / ^[x)(p{x)dx, (5.3.22) 

where (•, •) denotes the action of a tempered distribution on a test function (see 
Schwartz (1966)). Substituting the expression for (p given by (5.2.21) into the right 
hand side of (5.3.22) and using the fact that 

= V(0) = (2yr)- d / (5.3.23) 



gives 



(<p,4>) = -J ^(x)(l-G(t 1 / 2 x))t" 1 rfa(t)^ dx. 



We want to swap the order of integration here. This will be justified by Fubini's 
theorem if we can show that 

J \tp(x)\(l - G(t 1/2 x))t _1 da(t)^j dx < oo. (5.3.24) 

We defer the proof of (5.3.24) to Lemma 5.3.3 below and press on. Swapping the 
order of integration and recalling (5.3.23) yields 

(M) = -J q (/ ^(x)G(t 1/2 x) dx^j r 1 dait) 

= ~r ol md{ ^~ i/2) t_d/2_i da{t) 
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using Parseval's relation in the last line. Once again, we want to swap the order of 
integration and, as before, this is justified by Fubini's theorem if a certain integral 
is finite, specifically 

r ol im 1 d{it ~ i/2) ^) rd/2 ~ i Mt) < °°- (5 - 3 - 25) 

The proof of (5.3.25) will also be found in Lemma 5.3.3 below. After swapping the 
order of integration we have 

(M) = -f mm)d^ (5.3.26) 
which implies that <p = —H in lZ d \ {0}. 

Lemma 5.3.3. Inequalities (5.3.24) an d (5.3.25) are valid in Case I. 
Proof. For (5.3.24), we have 

J ^1 \^{x)\{l - G(* 1/2 x))* _1 da(*)^ dx 

< J (kJ^ \^(x)\\\x\\ 2 da(t)^j dx + J ^ ^(x)^- 1 da(t)j dx 

= K(a(l)-a(0)) J \^(x)\\\xf dx + ^ t^da^ (^J $(x)\dx^ 

< oo, 

recalling that ifj must enjoy faster than algebraic decay because ifj is a smooth 
function. 

For (5.3.25), the substitution r\ = £t -1 / 2 provides the integral 
I ■■= (J mvt 1/2 )\G(v) dn^j t- 1 da(t). 

Now there is a constant D such that \ip(y)\ < -D||y|| 2 for every y e 1Z d , because 
the support of ip is a closed subset of lZ d \ {0}. Hence 

I<J q D ^ G{rj) \\n\\ 2 dr^j da(t) + (2n) d G(0)U\\ oo r 1 da(t) 
< oo. 

The proof is complete. ■ 
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5.4. Lower bounds on eigenvalues 

Let (p:lZ d — > 1Z be a member of C and let (yj)j e z d be a zero-summing sequence. 
An immediate consequence of (5.2.26) is the equation 

E ViVwU-k) = (27r)- d I J £ fc-e 4 * |V(flde, (5.4.1) 

where = —H(£) for almost all £ G 7?. d and if is given by (5.2.27). Moreover, 
(5.2.6) is valid, that is 

E - *o = W / I E w e ' i€ 2(7 (o ^ ( 5 - 4 - 2 ) 

where a is given by (5.2.7). Applying (5.2.30), we have 

H0\= EN + 27rfc ) 



kez d 

r>00 



/>oo _ 

= / E G{t~ 1 / 2 {t + 2Tik)) t^ 2 - 1 da(t). 



(5.4.3) 



As in Section 2, we consider essential upper and lower bounds on a. Let us begin 
this study by fixing t > and analysing the function 



kez° 



By (5.3.14), we have 



where 



r 



^ l = l Ai (0) 



(5.4.4) 



(5.4.5) 



^■(x) = E A i ((x + 27rA;)t- 1 / 2 ), a; G TZ, j = l,...,d. (5.4.6) 

We now employ the following key lemma. 
Lemma 5.4.1. Let A G V and let 

E(x) = ^2k((x + 2ivk)t- 1/2 ), xeTZ. 
kez 

Then E is an even function and E(0) > E(x) > E(y) > E(ir) for every x and y 
in 1Z with < x < y < it. 
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Proof. The exponential decay of A and the absolute integrability of A imply that 
the Poisson summation formula is valid, which gives the relation 

E(x) = t 1 ' 2 A(H 1 / 2 )e ifcx , x G K. (5.4.7) 

kez 

Now the sequence a k := A(kt 1 / 2 ), k G Z, is an even, exponentially decaying 
Polya frequency sequence, that is every minor of the Toeplitz matrix (aj-k)j,ke.z 
is non-negative definite (and we see that this is a consequence of (5.3.7)). By 
a result of Edrei (1953), J2 keZ cikZ k is a meromorphic function on an annulus 
{z G C : 1/R < \z\ < R}, for some R > 1, and enjoys an infinite product expansion 
of the form 

^ GkZ = ° e 11 Z * °' (5 - 4 - 8) 

fce.z j=i v KJ /v ^ ; 

where C > 0, A > 0, < ctj, 8j < 1 and Xl^li a j + /^j < °°- Hence 

E(x)=Ct 1 / 2 e 2Xcosx Y\- — 3 - -f, (5.4.9) 

w 11 1 - 2B j cos a; + B 2 v ; 

Observe that each term in the product is an even function which is decreasing on 
[0, 2tt], which provides the required inequality. ■ 

In particular, Ej(x) > Ej(n) for j = 1, . . . , d, where Ej is given by (5.4.6). 

Hence 

r(£)>T(7re), £ g Tl\ (5.4.10) 

and applying (5.4.3) we get 

H0\ > W(7re)\, ieU d . (5.4.11) 
We now come to our principal result. 

Theorem 5.4.2. Let (yj)j e z d be a zero-summing sequence and let cp G C. Then 
we have the inequality 

VkyMj-k)\ > \<r(ve)\ ^ y 2 . (5.4.12) 

j,kez d jez d 
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Proof. Equation (5.4.2) and the Parseval relation provide the inequality 

£ VkVk<PU - k)\ > \a(ne)\(2n)- d j I £ y^tf d£ = \a(ne)\ £ y], 
j,kez d J[o,2-k]i jeZd jeZd 

as in inequality (5.2.10). ■ 

Of course, we are interested in showing that (5.4.12) cannot be improved, 
that is |cr(7re)| cannot be replaced by a larger number independent of (yj)j£z d - 
Recalling Proposition 5.2.2, this is true if a is continuous at ne. In fact, we can 
use Lemma 5.4.1 to prove that a is continuous everywhere in the set (0, 2it) d . We 
first collect some necessary preliminary results. 

Lemma 5.4.3. The function r given by (5.4-4) i> s continous for every t > and 
satisfies the inequality 

r(0 < r{rj) for < rj < £ < ne. (5.4.13) 

Furthermore, 

r ( 7re + f) = r (7re - for all £ E (-7T, n) d . (5.4.14) 
Proof. The definition of G, (5.4.5) and (5.4.7) provide the Fourier series 

r (£) = t d ' 2 Gikt 1 / 2 ^, £ e Tl d , (5.4.15) 

kez d 

and the exponential decay of G implies the uniform convergence of this series. 
Hence r is continuous, being the uniform limit of the finite sections of (5.4.15). 

Applying the product formula (5.4.5) and Lemma 5.4.1, we obtain (5.4.13) 
and (5.4.14). ■ 

Proposition 5.4.4. a is continuous on (0, 27r) d . 

2 

Proof. Equation (5.4.2) implies that Yljez d Vj e% ^ < 00 f° r almost every 

£ G [0,2-7r] d . Consequently, a is finite almost everywhere, by Lemma 5.2.1. Thus 
every non-empty open subset of [0, 27r] d contains a point at which a is finite. 
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Specifically, let 5 G (0,7r) and define the closed set Kg := [6,21V — 5] d . Thus the 
open set [0, 2ir} d \ Kg contains a point, r\ say, for which 

/•OO 

oo>\o-(r])\= V G((r t + 2nk)t- 1 / 2 )t- d / 2 - 1 da(t). (5.4.16) 
Jo kez* 

Let us show that a is continuous in Kg. To this end, choose any convergent se- 
quence (£ n )^Li in Kg and let £oo denote its limit. We must prove that 

lim <x(£ n ) = <r(foo)- 

n— »oo 

Now Lemma 5.4.3 and (5.4.16) supply the bound 

|(7(Cn)| < k(^)| < OO, 71=1,2,..., 

that is the functions 

{t ^ G((£n + 2nk)t- 1 / 2 )r d/2 - 1 da(t) : n = 1, 2, . . . } 

are absolutely integrable on [0, oo) with respect to the measure da. Moreover, 
they are dominated by the absolutely integrable function t h-> Ylkez d< ^^ r l + 
27vk)t~ 1 / 2 )t~ d / 2 ~ 1 . However, the continuity of r proved in Lemma 5.4.3 allows to 
deduce that 

lim £ G((£ n + 2^-^)^/2-1= £ G((^ oo + 2 7 rfc)t-V2 ) ^/2-i ? 

for all positive t. Thus the dominated convergence theorem implies that cr(£ n ) — > 
c(^oo) as n tends to infinity. Since 5 G (0, n) was arbitrary, we conclude that a is 
continuous in all of (0, 2n) d . ■ 

Corollary 5.4.5. Inequality (5.4-12) cannot be improved for cp G C if we can find 
a trigonometric polynomial P satisfying the conditions of Proposition 5.2.2 at the 
point ne. 

Proof. We simply apply Proposition 5.5.2. ■ 
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5.5. Total positivity and the Gaussian cardinal function 

This material is not directly related to the earlier sections of this chapter, but 
it does use a total positivity property to deduce an interesting fact concerning 
infinity norms of Gaussian distance matrices generated by infinite regular grids. 
Let A be a positive constant and let (p: 1Z — > 1Z be the Gaussian 

(p(x) = exp(— Xx 2 ), x G 1Z. (5.5.1) 

It is known (see Buhmann (1990)) that there exists a real sequence (ck)kez such 
that J2kez c 1 < 00 an d the function 7£ — )• 7£ given by 

X{x) = ^2c k <p(x - k), xeTZ, (5.5.2) 
kez 

satisfies the equation 

XU) = S 0j, j e Z. 

Thus x is the cardinal function of interpolation for the Gaussian radial basis 
function. 

Proposition 5.5.1. The coefficients (ck)kez of the cardinal function x alternate 
in sign, that is (—l) k Ck > for every integer k. 

Proof For each non-negative integer n, we let 

A n = (<p(j - k))l k= _ n . (5.5.3) 

Now A n is an invertible totally positive matrix, which implies that A~ x enjoys 
the "chequerboard" property, that is the elements of the inverse matrix satisfy 
(— iy +k (A~ 1 )jk > 0, for j, k = —n, . . . , n. In particular, if we let 

4 n) = (^Oofc, k = -n,...,n, (5.5.4) 
then (— l) k c^ > and the definition of A~ x provides the equations 

n 

c k fti - k ) = 5 oj, j = -n,...,n. (5.5.5) 

k=— n 
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In other words, the function \ n \ n^TZ defined by 



n 

» 



Xn(x) = c k ? f( x ~ xeTZ, (5.5.6) 



k= — n 



provides the cardinal function of interpolation for the finite set {— n, . . . , n}. 

Now Theorem 9 of Buhmann and Micchelli (1991) provides the following 
useful fact relating the coefficients of Xn and x : 



lim c[ = Ck, k E Z. 

Thus the property (— l) k c^ > implies the required condition (— l) k ct > 0. ■ 

We now consider the bi-infinite symmetric Toeplitz matrix A = (ip(j — 
k))j,kez as a bounded linear operator A: £ P (Z) — > £ P (Z) when p > 1. Thus A~ l = 
(cj-k)j,kez, where the (cj)j e z are given by (5.5.2), and a theorem of Buhmann 
(1990) provides the equation 

c k = (27T)" 1 / -L, e~ ik t df , k E Z, (5.5.7) 
Jo 

where 

Therefore, using standard results of Toeplitz operator theory (Grenander and 
Szego (1984)), we obtain the expression 

||A- 1 || 2 =max{^-:ee[0,27r]}. 

Applying Lemma 4.2.7, we get 

\\A-% = ^- = J2(- 1 ) kc K- ( 5 - 5 - 9 ) 

kTz 

But Proposition 5.5.1 and the symmetry of A provide the relations 

1 ii i = p-'iioo = \°*\ = E^- 1 )^' ( 5 - 5 - 10 )' 

kez kez 
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so that A~ x provides a nontrivial linear operator on l p (Z), for p = 1,2, and 
oo, whose norms agree on each of these sequence spaces. Further, we recall that 
log || A -1 ||p is a convex function of 1/p for p > 1, which is a consequence of the 
Riesz-Thorin theorem (Hardy et al (1952), pp. 214, 219). Hence we have proved 
the interesting fact that ||A _1 || P = ||A _1 ||i for all p > 1. 

In the multivariate case, the cardinal function is given by expressions anal- 
ogous to (5.5.2) and (5.5.7). Specifically, we let (f(x) = exp(— A||ir|| 2 ), x G 1Z d , and 
then x- 7t d ~^ Tt d is defined by 

X(x) = - AO, x G U d , (5.5.11) 

kez d 

where 

ci d) = (2n)- d [ -J—e- ik Ut, k = (h,...,k d )eZ d , (5.5.12) 



and 



^ (d) (0= E <£(£ + 27rA0. (5.5.13) 

kez d 

The key point is that <p is a tensor product of univariate functions, which implies 
the relation 

d 

0-^(0 = J] <r(&), £ = (6, ...,&) G (5.5.14) 
i=i 

where a is given by (5.5.8). Consequently the coefficients of the multivariate cardi- 
nal function are related to those of the univariate cardinal function by the formula 

d 

c k } = II % k = ••-Me (5.5.15) 

7=1 

In particular, the following corollary is an immediate consequence of Proposition 
5.5.1. 

Corollary 5.5.2. (-l) fc i+'"+ fc d c r) > q f or every integer k = (fci, . . . , kd) G 
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6 : Norm Estimates and Preconditioned Conjugate Gradients 



6.1. Introduction 

Let n be a positive integer and let A n be the symmetric Toeplitz matrix given by 

An = (<p(j ~ *))?,*=-» , (6.1.1) 

where (p: 1Z — > 1Z is either a Gaussian ((p(x) = exp(— \x 2 ) for some positive con- 
stant A) or a multiquadric (<p(x) = (x 2 + c 2 ) 1 / 2 for some real constant c). In this 
section we construct efficient preconditioners for the conjugate gradient solution 
of the linear system 

A n x = f, feU 2n +\ (6.1.2) 
when (p is a Gaussian, or the augmented linear system 

A n x + ey = f, 

(6.1.3) 

e T x = 0, 

when (p is a multiquadric. Here e = [1,1,..., 1] T G 7?. 2n+1 and y E 71. Sec- 
tion 6.2 describes the construction for the Gaussian and Section 6.3 deals with 
the multiquadric. Of course, we exploit the Toeplitz structure of A n to perform 
a matrix- vector multiplication in C(nlogn) operations whilst storing 0(n) real 
numbers. Further, we shall see numerically that the number of iterations required 
to achieve a solution of (6.1.2) or (6.1.3) to within a given tolerance is independent 
of n. 

Our method applies to many other radial basis functions, such as the inverse 
multiquadric (<p(x) = (x 2 + c 2 ) -1 / 2 ) and the thin plate spline (<p(x) = x 2 log \x\). 
However, we concentrate on the Gaussian and the multiquadric because they ex- 
hibit most of the important features of our approach in a concrete setting. Similarly 
we only touch briefly on the <i-dimensional analogue of (6.1.1), that is 



A ( n d) =(<Pti-k)) jtke[ . ntn] * 
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We shall still call A^f* a Toeplitz matrix. Moreover the matrix- vector multiplication 



. fc£[— n,n] d 



jE[-n,n 



d 



where || • || is the Euclidean norm and x = (xj) je [_ n n ]d, can still be calculated in 
0(N log N) operations, where N = (2n + l) d , whilst requiring O(N) real numbers 
to be stored. This trick is a simple extension of the Toeplitz matrix-vector multi- 
plication method when d = 1, but seems to be less familiar for d greater than one. 
This will be dealt with in detail in Baxter (1992c). 

6.2. The Gaussian 

Our treatment of the preconditioned conjugate gradient (PCG) method follows 
Section 10.3 of Golub and Van Loan (1989), and we begin with a general descrip- 
tion. We let n be a positive integer and A G TZ nXn be a symmetric positive definite 
matrix. For any nonsingular symmetric matrix P G TZ nXn and b G lZ n we can use 
the following iteration to solve the linear system PAPx = Pb. 

Algorithm 6.2.1. Choose any xq in 1Z n . Set tq = Pb — PAPxq and do = tq. 

For k = 0, 1, 2, ... do begin 
a k = rlr k /dlPAPd k 
Xk+i = Xk + a k dk 
r k+1 =r k - a k PAPd k 
b k = rl +1 r k+1 /rlr k 

d k+ i = r k+ i + b k d k 
Stop if ||r fc+1 || or Hdjt+ill is sufficiently small, 
end. 

In order to simplify Algorithm 6.2.1 define 

C = P 2 , Ck = Px k , r k = Pp k and S k = Pd k . (6.2.1) 
Substituting in Algorithm 6.2.1 we obtain the following method. 
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Algorithm 6.2.2. Choose any £o in 7t n ■ Set po = b — A^ , So = Cpo. 

For k = 0, 1, 2, . . . do begin 

a k = p\Cp k lb\Ab k 

= ik + dkO~k 

Pk+i = Pk~ ctkASk 
b k = pl +1 Cp k+1 /plCp k 

dk+i = Cpk+i + bk$k 
Stop if ||pfc + i|| or ||<5fc+i|| is sufficiently small, 
end. 

It is Algorithm 6.2.2 that we shall consider as our PCG method in this 
section, and we shall call C the preconditioner. We see that the only restriction 
on C is that it must be a symmetric positive definite matrix, but we observe that 
the spectrum of CA should consist of a small number of clusters, preferably one 
cluster concentrated at one. At this point, we also mention that the condition 
number of CA is not a reliable guide to the efficacy of our preconditioner. For 
example, consider the two cases when (i) CA has only two different eigenvalues, 
say 1 and 100, 000, and (ii) when CA has eigenvalues uniformly distributed in 
the interval [1,100]. The former has the larger condition number but, in exact 
arithmetic, the answer will be achieved in two steps, whereas the number of steps 
can be as high as n in the latter case. Thus the term "preconditioner" is sometimes 
inappropriate, although its usage has become standard. 

We can shed no light on the problem of constructing preconditioners for 
the general case.Accordingly, we let A be the matrix A n of (6.1.1) and let (p(x) = 
exp(— x 2 ). Thus A n is positive definite and can be embedded in the bi-infinite 
symmetric Toeplitz matrix 

Aoo = (<pU-k)) jMZ . (6.2.2) 

The classical theory of Toeplitz operators (see, for instance, Grenander and Szego 
(1984)) and the work of Section 4 provide the relations 

Sp A n c Sp A^ = [<t(7t), <j(0)} C (0, oo), (6.2.3) 
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where a is the symbol function 

cr(0 = J2^ + 2nk ^ ^n. (6.2.4) 
kez 

Further, Theorem 9 of Buhmann and Micchelli (1991) allows us to conclude that, 
for any fixed integers j and k, we have 

lim (A" 1 ),, k = (A^)jk. (6.2.5) 

n— >-oo 

It was equations (6.2.3) and (6.2.5) which led us to investigate the possibility of 
using some of the elements of A~ x for a relatively small value of n to construct 
preconditioners for An, where N may be much larger than n. Specifically, let us 
choose integers < m < n and define the sequence 

c 3 = ( A n 1 )j0, j = -m, . . . , m. (6.2.6) 
We now let Cn be the (2iV + 1) x (2A^ + 1) banded symmetric Toeplitz matrix 

f Cq ... Cm \ 



C 



N 



(6.2.7) 



V Cm ... Cq J 

We claim that, for sufficiently large m and n, Cn provides an excellent pre- 
conditioner when A = An in Algorithm 6.2.2. Before discussing any theoreti- 
cal motivation for this choice of preconditioner, we present an example. We let 
n = 64, m = 9 and iV = 32, 768. Constructing A n and calculating the elements 
{(A~ 1 )j : j = 0, 1, . . . , m} we find that 

/ 1.4301 x 10° \ 
-5.9563 x 10" 1 
2.2265 x 10" 1 
-8.2083 x 10" 2 
3.0205 x 10" 2 
-1.1112 x 10" 2 
4.0880 x 10" 3 
-1.5039 x lO- 3 
5.5325 x 10" 4 
V -2.0353 x 10" 4 / 



fc \ 

Cl 

\cj 



(6.2.8) 
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FIGURE 6.1: The symbol function for C^. 
Now Cn can be embedded in the bi-infinite Toeplitz matrix defined by 



(r ^ _ j Cj- k , \j -k\<m, 
lOooJjfc ~\0, \j-k\>m, 

and the symbol for this operator is the trigonometric polynomial 



(6.2.9) 



'0.(0= E c ^ t en - ( 6 - 2 - 10 ) 

J=—m 

In Figure 6.1 we display a graph of oc^ for < £ < 2ir, and it is clearly a 
positive function. Thus the relations 



Sp C N CS P C 00 = {a Coo (£) : £ G [0, 2tt]} C (0, oo) 
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imply that Cat is positive definite. Hence it is suitable to use Cat as the precon- 
ditioner in Algorithm 6.2.2. Our aim in this example is to compare this choice 
of preconditioner with the use of the identity matrix as the preconditioner. To 
this end, we let the elements of the righthandside vector b of Algorithm 6.2.2 be 
random real numbers uniformly distributed in the interval [—1, 1]. Applying Al- 
gorithm 6.2.2 using the identity matrix as the preconditioner provides the results 
of Table 6.1. Table 6.2 contains the analogous results using (6.2.7) and (6.2.8). 
In both cases the iterations were stopped when the residual vector satisfied the 
bound ||rfc+i||/||6|| < 10~ 13 . The behaviour shown in the tables is typical; we find 
that the number of steps required is independent of N and b. 



Iteration 


Error 




1 


2.797904 x 10 1 




10 


1.214777 x 10- 


2 


20 


1.886333 x 10" 


•6 


30 


2.945903 x 10" 


10 


33 


2.144110 x 10- 


•11 


34 


8.935534 x 10" 


•12 



Table 6.1: No preconditioning 



Iteration 


Error 




1 


2.315776 x 10" 


-l 


2 


1.915017 x 10" 


-3 


3 


1.514617 x 10" 


-7 


4 


1.365228 x 10" 


-11 


5 


1.716123 x 10" 


-15 



Table 6.2: Using (6.2.7) and (6.2.8) as the preconditioner 

Why should (6.2.7) and (6.2.8) provide a good preconditioner? Let us con- 
sider the bi-infinite Toeplitz matrix C^A^. The spectrum of this operator is given 
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by 

Sp Oo^oo = {<r Coo (£MO ■ £ e [0,2tt]}, (6.2.12) 

where a is given by (6.2.4) and oc x by (6.2.10). Therefore in order to concentrate 
Sp CooAoo at unity we must have 

^(CMO-l, ee[0,27r]. (6.2.13) 

In other words, we want oc^ to be a trigonometric polynomial approximating the 
continuous function 1/cr. Now if the Fourier series of 1/er is given by 

then its Fourier coefficients ('jj)jez are the coefficients of the cardinal function x 
for the integer grid, that is 

x (x) = Y / lM x -i), x^K, (6.2.15) 
jez 

and 

x(k)=5 Qk , kEZ. (6.2.16) 

(See, for instance, Buhmann (1990).) Recalling (6.2.5), we deduce that one way 
to calculate approximate values of the coefficients ('Yj)je.z is to solve the linear 
system 

A n c (n) = e°, (6.2.17) 

where e° = (Sjo) 1 j=- n ^ T<r n . This observation is not new; indeed Buhmann 
and Powell (1990) used precisely this idea to calculate approximate values of the 
cardinal function x- We now set 

Cj=cf\ 0<j<m, (6.2.18) 

and we observe that the symbol function cr for the Gaussian is a theta function 
(see Section 4.2). Thus a is a positive continuous function whose Fourier series is 
absolutely convergent. Hence 1/cr is a positive continuous function and Wiener's 
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lemma (Rudin (1973)) implies the absolute convergence, and therefore the uniform 
convergence, of its Fourier series. We deduce that the symbol function ac x can 
be chosen to approximate 1/ct to within any required accuracy. More formally we 
have the 

Lemma 6.2.3. Given any e > 0, there are positive integers m and no such that 



E 



>) e «€ _ i 



<e, £e[0,27r], 



j = -m 



for every n > no, where = (c^ )™=_ n is given by (6.2.17) 



Proof. The uniform convergence of the Fourier series for a implies that we can 
choose m such that 



j=-m 

By (6.2.5), we can also choose no such that |7j — | < e for j = — m, . . . , m and 
n > no- Then we have 



-(0 E 



>) e i* _ i 



< 



3=—m 



j = — m j = — m 

< e[l + (2m + l)<x(0)], 
remembering from Chapter 4 that < o~(n) < cr(£) < cr(0). Since e is arbitrary 
the lemma is true. ■ 



6.3. The Multiquadric 

The multiquadric interpolation matrix 

A = L(\\ Xj -x k \\)Y 

\ / J,k=l 

where <p(r) = (r 2 + c 2 ) 1 / 2 and (xj)" =1 are points in 1Z d , is not positive definite. 
We recall from Chapter 2 that it is almost negative definite, that is for any real 
numbers (yj)? = i satisfying J2yj =0 we have 

n 

E - x k\\) < °- (6.3.1) 

3,k=l 
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Furthermore, inequality (6.3.1) is strict whenever n > 2 and the points {xj)™ =1 are 
all different, and we shall assume this for the rest of the section. In other words, 
A is negative definite on the subspace < e >^, where e = [1, 1, . . . , 1] T G lZ n . 

Of course we cannot apply Algorithms 6.2.1 and 6.2.2 in this case. However 
we can use the almost negative definiteness of A to solve a closely related linearly 
constrained quadratic programming problem: 

minimize -£ T A£ — b T £ 

2^ (6.3.2) 

subject to e T £ = 0, 

where b can be any element of 1Z n . It can be shown that the standard theory of 
Lagrange multipliers guarantees the existence of a unique pair of vectors G 1Z n 
and 7]* G 1Z m satisfying the equations 

AC + er]* = b 

(6.3.3) 

and e T C = 0, 

where rj* is the Lagrange multiplier vector for the constrained optimization prob- 
lem (6.3.2). We do not go into further detail on this point because the nonsingu- 
larity of the matrix 

o) < 6 - 3 ' 4 > 

is well-known (see, for instance, Powell (1990)). Instead we observe that one way 
to solve (6.3.3) is to apply the following modification of Algorithm 6.2.1 to (6.3.2). 

Algorithm 6.3.1. Let P be any symmetric n x n matrix such that kerP = (e). 

Set xo = 0, r = Pb — PAPxq, do = r . 
For k = 0, 1, 2, ... do begin 
a k = rlr k /dlPAPd k 
x k +i = x k + a k d k 
r k+1 =r k - a k PAPd k 
b k = rl +1 r k+1 /rlr k 
d k+ i = rfc+i + b k d k 
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Stop if ||rfc+i|| or ||<ifc+i|| is sufficiently small, 
end. 

We observe that Algorithm 6.3. f solves the linearly constrained optimiza- 
tion problem 

minimize -x T PAPx — b T Px 

2 (6.3.5) 

rp 

subject to e x = 0. 
Moreover, the following elementary lemma implies that the solutions £*of (6.3.3) 
and x* of (6.3.5) are related by the equations = Px* . 

Lemma 6.3.2. Let S be any symmetric n x n matrix and let K = kerS. The 

S : K 1 - — > K 1 - is a bisection. In other words, given any b G K 1 - there is precisely 
one a G K 1 - such that 

Sa = b. (6.3.6) 
Proof. For any n x n matrix M we have the equation 

TZ n = kerM©Im M T . 
Consequently the symmetric matrix S satisfies 

TT = ker S © Im S, 

whence Im S = K^. Hence for every b G K 1 - there exists a G 1Z n such that 
Sa = b. Now we can write a = a + (3, where a G K 1 - and (3 G K are uniquely 
determined by a. Thus Sa = Sa = 6, and (6.3.6) has a solution. If a' G K 1 - also 
satifies (6.3.6), then their difference a — a' lies in the intersection K D K 1 - = {0}, 
which settles the uniqeuness of a. | 

Setting P = S and K = (e) in Lemma 6.3.2 we deduce that there is exactly 
one x* G (e) -1 such that 

PAPx* = Pb, 

and PAP is negative definite when restricted to the subspace (e)^. Follwing the 
development of Section 6.2, we define 

C = P 2 , Zk = Px k , and S k = Pd k , (6.3.7) 
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as in equation (6.2.1). However we cannot define p& by (6.2.1) because P is singular. 
One solution, advocated by Dyn, Levin and Rippa (1986), is to use the recurrence 
for (pfc) embodied in Algorithm 6.2.1 without further ado. 
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Algorithm 6.3.3a. Choose any £o in (e)" 1 . Set po = b — A£o and 5o = Cpo. 

For k = 0, 1, 2, . . . do begin 

a k = plCp k /5^A5 k 

= 6c + a kO~k 
Pk+1 = Pk — UkASk 

b k = pl +1 Cp k+1 /plCp k 
Sk+i = Cpt+i + bkSk 
Stop if ||pfc+i|| or is sufficiently small, 

end. 

However this algorithm is unstable in finite precision arithmetic, as we shall 
see in our main example below. One modification that sucessfully avoids instability 
is to force the condition 

Pk e (e)^, (6.3.8) 

to hold for all k. Now Lemma 6.3.2 implies the existence of exactly one vector pk G 
(e) 1 - for which Ppk = r^. Therefore, defining Q to be the orthogonal projection 
onto (e)- 1 , that is Q : x i— > x — e(e T x) / (e T e) , we obtain 

Algorithm 6.3.3b. Choose any ^ in Set p = Q(b — A^ ), 5 = Cp . 

For k = 0, 1, 2, . . . do begin 

a k = p^Cpk/S^ASk 

= ik + CLkSk 

Pk+i = Q(Pk - ctkASk) 
b k = pl +1 Cp k+ i/plCpk 

Sk+i = Cpk+i + bkd~k 
Stop if ||pfc + i|| or ||<5fc+i|| is sufficiently small, 
end. 

We see that the only restriction on C is that it must be a non-negative 
definite symmetric matrix such that kerC = (e). It is easy to construct such a 
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matrix given a positive definite symmetric matrix D by adding a rank one matrix: 

c = J _ (gW (6 . 3 . 9) 
e 1 JJe 

The Cauchy-Schwarz inequality implies that x T Cx > with equality if and only 
if x G (e). Of course we do not need to form C explicitly, since C : x (->• Dx — 
(e T Dx/e T De)De. Before constructing D we consider the spectral properties of 
Aoc = ((p(j - k)) jykeZ in more detail. 

A minor modification to Proposition 5.2.2 yields the following useful result. 
We recall the definition of a zero-summing sequence from Definition 4.3.1 and that 
of the symbol function from (5.2.7). 

Proposition 6.3.4. For every r\ G (0, 2ir) we can find a set {(Vj )jez '■ n = 
1,2,.. .} of zero-summing sequences such that 



j,kez jez 



Proof. We adopt the proof technique of Proposition 5.2.2. For each positive integer 
n we define the trigonometric polynomial 

n-l 
k=0 

and we recall from Section 4.2 that 



= ^r^ = \L n (0\\ (6-3.11) 
nsin 4/2 

where K n is the nth degree Fejer kernel. We now choose (yj )je.z to be the Fourier 
coefficients of the trigonometric polynomial £ t-> L n (^ — rj) sin£/2, which implies 
the relation 



X>fv* = sin2 e/2^n(e-"), 



and we see that (yj )jg.z is a zero-summing sequence. By the Parseval relation 
we have 

J>j n) ] 2 = (27T)" 1 / sin 2 £/2 K n (e-r/)^ (6.3.12) 



jez 
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and the approximate identity property of the Fejer kernel (Zygmund (1988), p. 
86) implies that 

r-2-K 

sin 2 77/2 = lim (27T)" 1 / sin 2 £/2 K n (i - 77) d£ 



= lim 5>^ 2 

jez 



Further, because a is continuous on (0, 2ir) (see Section 4.4), we have 

/>2tT 

sin 2 77/2 (7(7/) = lim (27T)" 1 / sin 2 £/2 K n (£ - 7/)<r(0 d£ 



lim V y) n) y { k n) ip(j - k), 



the last line being a consequence of (4.3.6). | 

Thus we have shown that, just as in the classical theory of Toeplitz op- 
erators (Grenander and Szego (1984)), everything depends on the range of val- 
ues of the symbol function a. Because a inherits the double pole that (p enjoys 
at zero, we have a: (0, 2tv) h-> (cr(7r),oo). In Figure 6.2 we display the function 
[0,2tt] 9£h> 1/(7(0- 

Now let m be a positive integer and let (dj)J l = _ rn be an even sequence of 
real numbers. We define a bi-infinite banded symmetric Toeplitz matrix by 
the equations 

(Ax,)* = I Ji;l- m ' (6-3.13) 



0, otherwise 

Thus (DooAoo)^ = ^(j - fc) where ^0*0 = Y^L-m d w{ x - 0- Further 

E yjv^U-k) = (27T)- 1 / ^^|^(£M0^- (6.3.14) 

Now the function £ h-> od^ (£)c(£) is continuous for £ e (0, 27r), so the argument 
of Proposition 6.3.4 also shows that, for every 77 G (0, 27r), we can find a set 
{(yi )je.H : n = 1, 2, . . . } of zero-summing sequences such that 



lim £ VTWW-*) Yfr T = ^(lWl)' (6-3-15) 
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Figure 6.2: The reciprocal symbol function 1/cr for the multiquadric. 
A good preconditioner must ensure that {od x (0 a (0 '• £ e (0, 27r)} is a 
bounded set. Because of the form of o~d x we have the equation 

m 
j=-m 

Moreover, as in Section 6.2, we want the approximation 

^(CMO-l, ee(0,27r), (6.3.17) 

and we need ctd^ to be a no n- negative trigonometric polynomial which is positive 
almost everywhere, which ensures that every one of its principal minors is positive 
definite. 
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Recalling Theorem 9 of Buhmann and Micchelli (1991), we let 



» 



j = -m, ...,m, 



(6.3.18) 



and to subtract a multiple of the vector [1, . . . , 1] T G 7^. 2m + 1 from (c < J n ^)JL_ rn to 
form a new vector (dj) 1 J i = _ m satisfying Y^j=- m dj = 0. Recalling that ~ 7j 
for suitable m and n, where 



" _1 (o = £7^* 



(6.3.19) 



and Yljez lj = ® (si nce °~ inherits the double pole of (p at zero), we hope to achieve 
(6.3.17). Fortunately, in several cases, we find that od^ is negative on (0, 2ir), so 
that ctd^ needs no further modifications. Unfortunately we cannot explain this 
lucky fact at present, but perhaps one should not always look a mathematical gift 
horse in the mouth. Therefore let n = 64 and m = 9. Direct calculation yields 



and we then obtain 



/c \ 

Cl 



/c \ 

Cl 



( -6.8219 x 10° \ 
4.9588 x 10° 
-2.0852 x 10° 
7.2868 x 10" 1 
-2.5622 x 10" 1 
8.8267 x 10" 1 
-3.1071 x 10" 2 
1.0626 x 10- 2 
-3.7923 x 10" 3 
V 1.2636 x 10-3 / 



-6.8220 x 10° \ 

4.9587 x 10° 
-2.0852 x 10° 

7.2863 x 10" 1 
-2.5626 x 10" 1 

8.8224 x 10" 1 
-3.1113 x 10" 2 

1.0583 x 10" 2 
-3.8350 x 10" 3 

1.2210 x 10" 3 / 



(6.3.20) 



(6.3.21) 



110 



Preconditioned conjugate gradients 

Figures 6.3 and 6.4 display the functions od^ and £ >->■ <JD x {i)/ sin 2 (£/2) on the 
domain [0, 2n] respectively. The latter is clearly a positive function, which implies 
that the former is positive on the open interval (0, 27r). 
Thus, given 

A N =(<p(j-k)) 

\ / j,k=-N 

for any N > n, we let Dn be any (2N + 1) x (2N + 1) principal minor of and 
define the preconditioner CV by the equation 



(D N e)(D N 



C N = D N - V JV T n » ( 6 - 3 - 22 ) 

where e = [1, . . . , 1] T G TZ 2N+1 . We reiterate that we actually compute the matrix- 
vector product Cnx by the operations x >->■ Datx — {e T D / e T D ^ e)e rather than 
by storing the elements of CV in memory. 

Cn provides an excellent preconditioner. Tables 6.3 and 6.4 illustrate its 
use when Algorithm 6.3.3b is applied to the linear system 

A N x + ey = b, 



e T x = 0, 



(6.3.23) 



when N = 2, 048 and iV = 32, 768 respectively. Here y E TZ, e = [1, . . . , 1] T € 
1Z 2N+1 and b G 1Z 2N+1 consists of pseudo-random real numbers uniformly dis- 
tributed in the interval [—1,1]. Again, this behaviour is typical and all our numer- 
ical experiments indicate that the number of steps is independent of N. We remind 
the reader that the error shown is but that the iterations are stopped when 

either or ||5fe+i|| is l ess than 10 _13 ||6||, where we are using the notation of 

Algorithm 6.3.3b. 

It is interesting to compare Table 6.3 with Table 6.5. Here we have chosen 
ni= 1, and the preconditioner is essentially a multiple of the second divided 
difference preconditioner advocated by Dyn, Levin and Rippa (1986). Indeed, we 
find that do = 7.8538 and d\ = <i_i = —3.9269. We see that its behaviour is clearly 
inferior to the preconditioner generated by choosing m = 9. Furthermore, this is 
to be expected, because we are choosing a smaller finite section to approximate 
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the reciprocal of the symbol function. However, because (JD^iO is a multiple of 
sin 2 £/2, this preconditioner still possesses the property that {o'D oc (O a (0 '■ £ £ 
(0, 2tt)} is a bounded set of real numbers. 



Iteration Error 



1 


3 


.975553 


X 


10 4 




2 


8 


.703344 


X 


10" 


•l 


3 


2 


.463390 


X 


io- 


-2 


4 


8 


.741920 


X 


10" 


•3 


5 


3 


.650521 


X 


10" 


-4 


6 


5 


.029770 


X 


10" 


•6 


7 


1 


.204610 


X 


10" 


•5 


8 


1 


.141872 


X 


10" 


-7 


9 


1 


.872273 


X 


10" 


-9 


10 


1. 


.197310 


X 


10" 


•9 


11 


3 


.103685 


X 


10" 


-11 



Table 6.3: Preconditioned CG - m = 9, n = 64, iV = 2, 048 



Iteration 


Error 








1 


2.103778 


X 


10 5 




2 


4.287497 


X 


10° 




3 


5.163441 


X 


10" 


i 


4 


1.010665 


X 


10" 


i 


5 


1.845113 


X 


io- 


3 


6 


3.404016 


X 


io- 


3 


7 


3.341912 


X 


io- 


5 


8 


6.523212 


X 


10" 


7 


9 


1.677274 


X 


10" 


5 


10 


1.035225 


X 


io- 


8 


11 


1.900395 


X 


io- 


10 



Table 6.4: Preconditioned CG - m = 9, n = 64, = 32, 768 
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It is also interesting to compare the spectra of C n A n for n = 64 and m = 1 
and m = 9. Accordingly, Figures 6.5 and 6.6 display all but the largest nonzero 
eigenvalues of C n A n for m = 1 and m = 6 respectively. The largest eigenvalues 
are 502.6097. and 288.1872, respectively, and these were omitted from the plots in 
order to reveal detail at smaller scales. We see that the clustering of the spectrum 
when m = 9 is excellent. 

Iteration Error 



1 


2, 


645008 


X 


10 4 




10 


8. 


.632419 


X 


10° 




20 


9. 


.210298 


X 


10" 


i 


30 


7. 


.695337 


X 


10" 


i 


40 


3 


.187051 


X 


10" 


5 


50 


5. 


.061053 


X 


io- 


7 


60 


7. 


.596739 


X 


10" 


9 


70 


1. 


.200700 


X 


10" 


10 


73 


3 


,539988 


X 


10" 


11 


74 


1. 


,992376 


X 


10" 


11 



Table 6.5: Preconditioned CG - m = 1, n = 64, N = 8, 192 

The final topic in this section demonstrates the instability of Algorithm 
6.3.3a when compared with Algorithm 6.3.3b. We refer the reader to Table 6.6, 
where we have chosen m = 9, n = N = 64, and setting b = [1, 4, 9, ... , N 2 ] T . The 
iterations for Algorithm 6.3.3b, displayed in Table 6.7, were stopped at iteration 
108. For Algorithm 6.3.3a, iterations were stopped when either ||pfc+i|| or ||5fc+i|| 
became smaller than 10 _13 ||6||. It is useful to display the norm of \\8k\\ rather than 
\\pk\\ in this case. We see that the two algorithms almost agree on the early iter- 
ations, but that Algorithm 6.3.3a soon begins cycling, and no convergence seems 
to occur. Thus when can leave the required subspace due to finite precision 
arithmetic, it is possible to attain non-descent directions. 
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Iteration 


\\5k - 6.3.3a 


\\5k - 6.3.3b 


1 


4.436896 x 


10 4 


4.436896 x 10 4 


2 


2.083079 x 


10 2 


2.083079 x 10 2 


3 


2.339595 x 


10° 


2.339595 x 10° 


4 


1.206045 x 


10" 1 


1.206041 x 10" 1 


5 


1.698965 x 


10" 3 


1.597317 x 10" 3 


6 


6.537466 x 


10" 2 


6.512586 x 10" 2 


7 


1.879294 x 


10" 4 


9.254943 x 10" 6 


8 


2.767714 x 


10" 2 


1.984033 x 10" 7 


9 


3.453789 x 


10" 4 




10 


1.914126 x 


lO- 3 




20 


4.628447 x 


10" 1 




30 


3.696474 x 


io-° 




40 


8.061922 x 


10+ 3 




50 


2.155310 x 


10° 




100 


3.374467 x 


10" 1 




101 


1.121903 x 


10° 




102 


1.920517 x 


10" 1 




103 


3.772007 x 


10" 2 




104 


3.170231 x 


10" 2 




105 


2.612073 x 


10" 1 




106 


2.236274 x 


10° 




107 


8.875137 x 


lO" 1 




108 


1.823607 x 


10" 1 





Table 6.5: Algorithms 6.3.3a & b - m = 1, n = 64, N = 64, b = [1, 4, . . . , N 2 ] T ■ 
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Preconditioned conjugate gradients 




0.99 



Figure 6.6: The spectrum of C n A n for m = 9 and n = 64. 
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7 : On the asymptotic cardinal function for the multiquadric 



7.1. Introduction 

The radial basis function approach to interpolating a function /: 1Z d — > 1Z on the 
integer lattice Z d is as follows. Given a continuous univariate function (p: [0, oo) — > 
1Z, we seek a cardinal function 

X (x)=J2 a M\\x-j\\), x£K d : (7.1.1) 



that satisfies 



Therefore 



x (k) = 6 0tk , k e Z d . 



if(x)= fU)x(x-j), xen d , (7.1.2) 

jez d 

is an interpolant to / on the integer lattice whenever (7.1.2) is well defined. Here 
|| • || is the Euclidean norm on lZ d . This approach provides a useful and flexible 
family of approximants for many choices of (p, but here we concentrate on the 
Hardy multiquadric <p c (r) = (r 2 + c 2 ) 1 / 2 . For this function, Buhmann (1990) has 
shown that a cardinal function Xc exists and its Fourier tranform is given by the 
equation 

* c (0 = ^ Wtto hiv * G ^' (7 - L3) 

l^kez* VM||4 + 27r/c||) 
where {<£ c (||£||) : £ e ^ s the generalized Fourier transform of {(/? c (||x||) : 

x G 7Z d }. Further, Xc possesses a classical Fourier transform (see Jones (1982) or 
Schwartz (1966)). In this chapter, we prove that Xc enjoys the following property: 



which sheds new light on the approximation properties of the multiquadric asc-> 
oo. For example, in the case d = 1, (7.1.4) implies that lim^oo Xc(x) = sinc(x), 
providing a perhaps unexpected link with sampling theory and the classical theory 
of the Whittaker cardinal spline. Further, our work has links with the error analysis 
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of Buhmann and Dyn (1991) and illuminates the explicit calculation of Section 4 
of Powell (1991). It may also be compared with the results of Madych and Nelson 
(1990) and Madych (1990), because these papers present analogous results for 
polyharmonic cardinal splines. 

7.2. Some properties of the multiquadric 

The generalized Fourier transform of <p c is given by 

MUW) = -K-\2Kc/\\m {d+l)/2 K {d+lV2 {c\\m, (7-2.1) 

for nonzero £ G lZ d (see Jones (1982)). Here {K„(r) : r > 0} are the modified 
Bessel functions, which are positive and smooth in 1Z + , have a pole at the origin, 
and decay exponentially (see Abramowitz and Stegun (1970)). There is an integral 
representation for these modified Bessel functions (Abramowitz and Stegun (1970), 
equation 9.6.23) which transforms (7.2.1) into a highly useful formula for (p c : 

/oo 
eM-cxU\\)(x 2 ~ l) d/2 dx, (7.2.2) 

where Ad = 7r d / 2 /r(l + d/2). A simple consequence of (7.2.2) is the following 
lemma, which bounds the exponential decay of (p c . 

Lemma 7.2.1. // ||f || > \r)\ > 0, then 

i^(iieii)i<exp[- c (iieii-iHi)]i^ c (iHi)|. 

Proof. Applying (7.2.2), we obtain 

/oo 
eM~cx(U\\ - IMI)] exp(-c:r|M|) (x 2 - l) d / 2 dx 

<eM-c(M\-H\))m\\r,\\)l 
providing the desired bound. | 

We now prove our main result. We let I:TZ — > 1Z be the characteristic 
function of the cube [— tt, n] d , that is 

lu to, ii[-^M d - 

120 



On the asymptotic cardinal function for the multiquadric 
Proposition 7.2.2. Let £ be any fixed point of1Z d . We have 



lim *c(£) = /(£), 

c— >oo 

if llClloo 7^ ^ ^<*£ 2s £ ^oes He ^ e boundary of [— 7r,7r] d . 

Proof. First, suppose that £ ^ [— 7r,7r] d . Then there exists a nonzero integer /c 
such that ||£ + 27r/c || < ||£||, an d Lemma 7.2.1 provides the bounds 

l^cdlClDI < exp[-c(||e|| - ||e + 2nk G \\)]\M\t + ^k \\)\ 

< eX p[-c(||£|Ml£ + 27rM)] E + 27rfc||)|. 

kez d 

Thus, applying (7.1.3) and remembering that ip c does not change sign, we have 



< xc(0 < exp[-c(||e|| - ||e + 2^11)], e £ [-tt, 7r] d . 



(7.2.3) 



The upper bound of (7.2.3) converges to zero as c — > oo, which completes the proof 
for this range of £. 

Suppose now that £ G (— 7r, 7r) d . Further, we shall assume that £ is nonzero, 
because we know that % c (0) = 1 for all values of c. Then ||£ + 2irk\\ > ||£||, for 
every nonzero integer k G Z d . Now (7.1.3) provides the expression 



* c (o = (i+ E 

kez d \{o} 



¥c(U + 27Tk\ 



0c(U\\) 



)■■ 



(7.2.4) 



We shall show that 



lim > 



fc£.E d \{0} 



^c(||£ + 27T/c| 



MUW) 



= 0, £G(-7r,7r) d , 



which, together with (7.2.4), implies that lim^oo Xc(0 = 1- 
Now Lemma 7.2.1 implies that 



fce.z d \{o} 



<p c (U + 2irk\ 



<Pc(U\\) 



< E exp[-c(||£ + 27rA;| 
kez d \{o} 



(7.2.5) 



IICII)], (7-2.6) 
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and each term of the series on the right converges to zero as c — > oo, since ||£ + 
2nk\\ > ||£|| for every nonzero integer k. Therefore we need only deal with the tail 
of the series. Specifically, we derive the equation 

hm V exp[-c(||e + 27rA;||-||e||)] = 0, (7.2.7) 

c— >oo — ' 

ll*ll>2||e|| 

where e = [1, 1, . . . , 1] T . Now, if \\k\\ > 2||e||, then 

||e + 27rfc||-||e||>27r(||fc||-||e||)>7r||fc||, 
remembering that we have ||£|| < 7r||e||. Hence 

exp[-c(||e + 27rA;||-||e||)]< V exp(-7rc||fc||). (7.2.8) 



^ - -M'i • • is ■ — m ism ' j - 

|fc||>2||e|| Hfc||>2||e| 



It is a simple exercise to prove that the series S||fe||>2||e|| ex P( — ^ll^ll) ^ s convergent. 
Therefore, given any e > 0, there exists a positive number R > 1 such that 

exp(-7r||fc||) < e. 

||fc||>2A||e|| 

Consequently, when c > \R] we have the inequality 

J2 exp(-7rc||fc||) < exp(-7r||A:||) < e, 

||fe||>2||e|| l|fc||>2fl||e|| 

which establishes (7.2.5). The proof is complete. | 



7.3. Multiquadrics and entire functions of exponential type 7r 

Definition 7.3.1 Let / e L 2 {TZ d ). We shall say that / is a. function of exponential 
type A if its Fourier transform / is supported by the cube [—A, A] d . We shall denote 
the set of all functions of exponential type A by EA{7i d )- 

We remark that the Paley- Wiener theorem implies that / may be extended 
to an entire function on C d satisfying a certain growth condition at infinity (see 
Stein and Weiss (1971), pages 108ff), although we do not need this result. 
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Lemma 7.3.2. Let f G E 7r (Tl d )r\L 2 (Tl d ) be a continuous function. Then we have 
the equation 

E /(* + 2wk) = /(*) «*pH*0. (7-3.1) 

kez d kez d 

the second series being convergent in L 2 (lZ d ). 
Proof. Let 

9(0= E ht+ 2 **>)> ^ nd - 

kez d 

At any point £ G 1Z d , this series contains at most one nonzero term, because of 
the condition on the support of /. Hence g is well defined. Further, we have the 
relations 

< oo, 



/ l<7(£)| 2 ^= / |/(OI 2 d6 

J\—K,Tr] d Jn d 



since the Parseval theorem implies that / is an element of L 2 (1Z d ). Thus g G 
L 2 ([— 7r, 7r] d ) and its Fourier series 

= E 9kexp(ik^), 
kez d 

is convergent in L 2 ([— n, tt] rf ) . The Fourier coefficients are given by the expressions 
g k = (2n)~ d [ /(£) exp(-ifcO rf£ = (2vr)- d / /(£) exp(-tfcO d£ = f(-k), 

J[--K,TT] d Jn d 

where the final equation uses the Fourier inversion theorem for L 2 (1Z d ) . The proof 
is complete. | 

We observe that an immediate consequence of the lemma is the convergence 
of the series Ylkez d t/(^)] 2 ' by the Parseval theorem. 

For the following results, we shall need the fact that Xc G L 2 {1Z d ), which is 
a consequence of the analysis of Buhmann (1990). 

Lemma 7.3.3. Let f G E^TZ' 1 ) D L 2 (1Z d ) be a continuous function. For each 
positive integer n, we define the function 

SffiO = I E /(*) exp(-z^) I xc(£), £ G U d . (7.3.2) 

\||fc||i<n / 

Then {S™f : n = 1, 2, . . .} forms a Cauchy sequence in L 2 (1Z d ). 
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Proof. Let Q n : 1Z d — >■ 1Z be the trigonometric polynomial 

Q»(0= E /(*)exp(-ifcO, (7-3.3) 

||fc||i<n 

so that S™f(£) = Qn(C)Xc(C)- 1^ is a consequence of Lemma 7.3.2 that this se- 
quence of functions forms a Cauchy sequence in L 2 ([—n,7c] d ). Indeed, we shall 
prove that for m > n we have 

\\SFf - SjTfWw) < \\Qm - Qn\\ L *([-«,«]*), (7-3.4) 

so that the sequence of functions {S™f : n = 1, 2, . . .} is a Cauchy sequence in 

L 2 (n d ). 

Now Fubini's theorem provides the relation 

\\sFf - sjTfWhv*) = J nd \Qm(0 - Qn(d)\ 2 *&) # 



•/[-Tr.Tr]" ^ 6 ^ d y 



(7.3.5) 

However, (7.1.3) gives the bound 

xl{i + 2*i) = E <p 2 m + 2 ^ii)/( E + 27r *H)) 2 

ze.z d ze.H d fc6^ d (7.3.6) 
<1, 

which, together with (7.3.5), yields inequality (7.3.4). | 
Thus we may define 

= Xc(0 E /(*) «p(-ifcO, (7-3.7) 

fce.E d 

and the series is convergent in L 2 (1Z d ). Applying the inverse Fourier transform 
term by term, we obtain the useful equation 

SJ(x)= E f{k)Xc{x-k), xell d . 
kez d 

Theorem 7.3.4. Let f G E n (TZ d ) n L 2 {1Z d ) be a continuous function. We have 

lim S c f(x) = f(x), 

c— >oo 

and the convergence is uniform on TZ d . 
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Proof. We have the equation 

SJ(x) - f(x) = (2n)- d [ £ /(* + 2nk) ( Xc (0 ~ I(£))exp(ix£) df. 

Thus we deduce the bound 
|5 c /(x)-/(x)| 



<(2n)~ d [ £ |xc(^ + 27rfc)-/(^ + 27rfc) 

= (27r)- d / i/(oi ( i - xc(o + E + 2yrfc ) ] 

(7.3.8) 

using the fact that Xc is non-negative, and we observe that this upper bound is 
independent of x. Therefore we prove that the upper bound converges to zero as 
c — > oo. 

Applying (7.1.3), we obtain the relation 

£ Xc(£ + 27rfc) = l-xc(0, (7-3-9) 
fces d \{o} 



whence 



|5 c /(x) - /(z)| < 2(2yr)- d / 1/(01(1 - Xc(0) (7.3.10) 

J[-7T,7r] d 

Now / G L 2 ([— 7T, 7r] d ) implies / G -^ 1 ([— 7r, 7r] ), by the Cauchy-Schwartz inequal- 
ity. Further, Proposition 7.2.2 gives the limit lim^oo Xc(Q = 1, for £ G (— 7r,7r) d , 
and we have < 1 — % c (0 < 1, by (7.1.3). Therefore the dominated convergence 
theorem implies that 

hm(27r)- d / |/(0|(1 -Xc(O)df = 0. 
The proof is complete. | 
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On the asymptotic cardinal function for the multiquadric 
7.4. Discussion 

Section 4 of Powell (1991) provides an explicit calculation that is analogous to the 
proof of Theorem 7.3.4 when f(x) = x 2 . Of course, this function does not satisfy 
the conditions of Theorem 7.3.4. Therefore extensions of this result are necessary, 
but the final form of the theorem is not clear at present. 

Theorem 7.3.4 encourages the use of large c for certain functions. Indeed, 
it suggests that large c will provide high accuracy interpolants for univariate func- 
tions that are well approximated by integer translates of the sine function. Thus, in 
exact arithmetic, a large value of c should be useful whenever the function is well 
approximated by the Whittaker cardinal series. However, we recall that the linear 
systems arising when c is large can be rather ill-conditioned. Indeed, in Chapter 
4 we proved that the smallest eigenvalue of the interpolation matrix generated by 
a finite regular grid converges to zero exponentially quickly as c — > oo. We refer 
the reader to Table 4.1 for further information. Therefore special techniques are 
required for the effective use of large c. 
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Conclusions 
8 : Conclusions 

There seems to be no interest in using non-Euclidean norms for radial basis 
functions at present, possibly because of the poor approximation properties of 
the £ 1 -norm || • ||i reported by several workers. Thus Chapter 2 does not seem 
to have any practical applications yet. However, it may be useful to use p-norms 
(1 < p < 2), or functions of p- norms, when there is a known preferred direction 
in the underlying function, because radial basis functions based on the Euclidean 
norm can perform poorly in this context. On a purely theoretical note, we observe 
that the construction of Section 2.4 can be applied to any norm enjoying the 
symmetries of the cube. 

The greatest weakness - and the greatest strength - of the norm estimates 
of Chapters 3-6 lies in their dependence on regular grids. However, we note that 
the upper bounds on norms of inverses apply to sets of centres which can be 
arbitrary subsets of a regular grid. In other words, contiguous subsets of grids are 
not required. Furthermore, we conjecture that a useful upper bound on the norm 
of the inverse generated by an arbitrary set of centres with minimal separation 
distance 5 (that IS 1 1 X j 

| > 5 > if j ^ k) will be provided by the upper 
bound for the inverse generated by a regular grid of spacing 8. 

Probably the most important practical finding of this dissertation is that 
the number of steps required by the conjugate gradient algorithm can be indepen- 
dent of the number of centres for suitable preconditioners. We hope to discover 
preconditioners with this property for arbitrary sets of centres. 

The choice of constant in the multiquadric is still being investigated (see, for 
instance, Kansa and Carlson (1992)). Because the approximation of band-limited 
functions is of some practical importance, our findings may be highly useful. In 
short, we suggest using as large a value of the constant as the condition number 
allows. Hence there is some irony in our earlier discovery that the condition number 
of the interpolation matrix can increase exponentially quickly as the constant 
increases. 

Let us conclude with the remark that radial basis functions are extremely 
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rich mathematical objects, and there is much left to be discovered. It is our hope 
that the strands of research initiated in this thesis will enable some of these future 
discoveries. 
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